Non-Canonical NF-κB Signaling Stratifies LGG into Subtypes with Distinct Molecular and Cellular Characteristic and Survival Expectancy

Introduction

Low-grade glioma (LGG) accounts for approximately 25% of primary brain tumors. Despite advances in treatment development, more than half of LGG tumors progress into therapy-resistant high-grade gliomas with poor prognosis.¹ Several studies have revealed that LGG is genetically heterogeneous.² Histopathological classification combined with the analysis of recurrently mutated genes such as isocitrate dehydrogenase 1 (IDH1), IDH2, TP53, EGFR, and ATRX is used to determine the prognosis of LGG patients.^3,4 In addition, co-deletion of chromosome arms 1p and 19q (1p/19q codeletion) and O-6-methylguanine-DNA methyltransferase (MGMT) promoter methylation status have been identified as prognostic factors for LGG.⁵ However, LGG patients with the same risk factors have different prognosis based on these markers, indicating that more prognostic factors remain unrevealed.

The tumor microenvironment (TME) mediates cancer progression, responses to therapy and prognosis. Tumor-infiltrating immune cells (TIICs) are critical components of the TME and are associated with clinical outcomes.⁶ TIICs mainly include tumor-infiltrating lymphocytes (TILs), tumor-associated macrophages (TAMs), natural killer (NK) cells, dendritic cells, and myeloid-derived suppressor cells (MDSCs). TAMs consist of a small portion of M1 macrophages and a large portion of M2 macrophages. Interferon-γ (IFN–γ) and/or lipopolysaccharides (LPS) induce M1 macrophage polarization, while IL-4, IL-13 or TGFb1 induces M2 macrophage polarization.⁷ M1 macrophages promote the generation of T helper type 1 lymphocytes, thus fighting tumors, while M2 macrophages promote tumor progression. M2 macrophages are associated with poor prognosis in high-grade glioma.⁸ TILs are usually associated with favorable outcomes; in contrast, immunosuppressive cells such as regulatory T (Treg) cells, MDSCs and M2 macrophages are associated with poor outcomes.⁹ These findings suggest that TIICs might be the determining factors for LGG heterogeneity. Therefore, elucidating the underlying TIIC mechanisms is essential to provide actionable targets for immunotherapy to achieve precise treatment for LGG.

The NF-κB signaling is a major pathway that mediates TIICs.¹⁰ The NF-κB family has five core members: RelA (p65), RelB, Rel, NF-κB1 (p50 and its precursor p105), and NF-κB2 (p52 and its precursor p100). These core members form dimers to regulate two distinct NF-κB signaling activation mechanisms: the canonical NF-κB signaling (p50/RelA) and the non-canonical NF-κB signaling (p52/RelB). An atypical NF-κB signaling lacks both RelA and RelB. Inhibiting NF-κB activation in TAMs converts them from M2 to M1 macrophages.¹¹ IL-1β activates MDSCs via NF-κB signaling.¹² NF-κB activation is strongly associated with T cell infiltration in human lung cancer.¹³ Although less studied comparing with canonical NF-κB signaling, the non-canonical NF-κB signaling has been known to regulate immune response as well as the development of immune related organs. For instance, non-canonical NF-κB involves thymus development by regulating the development of thymic epithelial cells.¹⁴ Non-canonical NF-κB also plays a role in the development of secondary lymphoid organs, including spleen and lymph node.¹⁵ In terms of immune response, non-canonical NF-κB signaling participates in the generation as well as maintenance of effector and memory T cells.^16,17 Dendritic cells (DCs) is the primary antigen presenting cells linking innate immunity and adaptive immunity. Its maturation has been also found to be regulated by non-canonical NF-κB signaling.¹⁸ While the role of non-canonical NF-κB signaling in regulating immunity is being uncovered, the interaction between dysregulated non-canonical NF-κB signaling and immune environment in cancerous context remains unknown.

Here, we investigated the activity of non-canonical NF-κB signaling in LGG by leveraging multi-omics data from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA).¹⁹ We found that non-canonical NF-κB signaling can divided LGG into molecularly distinct groups. Particularly, non-canonical NF-κB-high group enriched with LGG featured by oncogenic mutations and immunosuppressive environment, consistently aligned with the poor survival of LGG patients in this groups. Our analysis provided potential patient stratification strategy for better therapeutic strategy as well as prognosis.

Materials and Methods

Data Retrieval and Preprocessing

The RNA-Seq data, somatic mutation data and clinical information of LGG tissue samples were downloaded from The Cancer Genome Atlas (TCGA) and the Chinese Glioma Genome Atlas (CGGA) data portal. Samples that met the following criteria were included in our study: 1) primary tumor samples; 2) both somatic mutation and expression data were available; 3) patients’ clinical information included both prognosis and previous biomarker test results is available. The expression data were normalized as fragments per kilobase million (FPKM). Somatic mutation data was downloaded in MAF format. The clinical data were integrated and downloaded through cBioPortal (http://www.cbioportal.org).

Estimation of NF-κB Signaling Activity

The gene sets of canonical NF-κB signaling, non-canonical NF-κB signaling and atypical NF-κB signaling were retrieved from the Molecular Signatures Database (MSigDB).²⁰ The pathway enrichment score was calculated using single sample gene set enrichment analysis (ssGSEA, v4).²⁰ Survival analysis was conducted using Kaplan-Meier estimator. Patients stratification was done using X-Tile.²¹

Characterization of Tumor-Infiltrating Immune Cells (TIICs) and Immunogenicity Estimation

ESTIMATE (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data: https://bioinformatics.mdanderson.org/estimate/ https://www.gsea-msigdb.org/gsea/index.jsp) was used to estimate the relative fraction of tumor cells, stromal cells and immune cells in LGG sample. ImmunecellAI (http://bioinfo.life.hust.edu.cn/ImmuCellAI) was used for immune cell decomposition.

Tumor mutation burden (TMB) is defined as the average amount of somatic mutations per Mb within the protein coding region. Somatic mutations were counted only if they meet the following criteria: 1) mutant allele frequency > 5%; 2) mutant does not included in common_dbSNP of the dbSNP knowledgebase; 3) population frequency of mutant < 1% in 1000 Genomes Project. Immunophenogram^22,23 was used to analyze the mRNA expression data of the TCGA-LGG dataset to estimate Immunophenotypescore (IPS), Briefly, immunogenicity considers four categories of immune-relevant gene sets representing effector cells (ECs), suppressive cells (SCs), major histocompatibility complex molecules (MHCs), and immunomodulators or checkpoints (CPs), The IPS was scaled to 0–10 based on weighted averaged Z-scores of the determinants from the above four categories and increased with higher immunogenicity, as previously described.²² Similarly, T cell–inflamed GEP was used to estimate the immunogenicity of LGG samples.²⁴

Survival Analysis and Construction of Prognostic Model

Survival analysis was performed with the R package survival analysis. For univariate survival analysis, Log rank test was performed to determine statistical difference. In terms of multivariate survival analysis, Cox proportional hazards model was applied to determine the independence of prognostic predictor by adjusting potential confounders. P value less than 0.05 was considered as statistically significant.

To constructing the prognostic model, TCGA-LGG patients were randomly divided into a training set (n=273) and an internal testing set (n=137) at ratio of 2:1. Through univariate regression, 91 genes (Supplementary Table 2) were prioritized from 165 genes (Supplementary Table 1) that were initially collected from MSigDB. To better generalize our prognostic model, we performed least absolute shrinkage and selection operator (LASSO) regression to select 9 genes (Supplementary Table 3) with greater prognostic power. To account for potential confounders, such as age, sex, grade, MGMT promoter methylation status, 1p19q status, and IDH, we further used multivariate Cox regression to shrink the final gene list to nine genes whose expression levels were used to construct a risk score system based on the following formula:

where “n” is the total number of prognostic genes; β_i and E_i represent the LASSO regression coefficient of gene i and log2 transformed relative expression value of gene i. Based on the median risk score, LGG patients were divided into a high-risk group and low-risk group. A multivariate Cox model was used to assess the impact of the model on the survival of LGG patients. To validate the prognostic capability of our model, an independent dataset was retrieved from the CGGA dataset and subjected to our model after preprocessing. The ROC was generated using R package “timeROC”.

Statistical Analysis

R language (v3.6.1) was used for statistical analysis. Between groups comparison was conducted using either Mann–Whitney U-test or Kruskal–Wallis test depending on the number of groups. Correlation between continuous variables was measured using Spearman correlation coefficient. Considering the large sample size subjected to statistical analysis, in addition to P-value (*p<0.05, **p<0.01 and ***p<0.001), fold change or Pearson coefficient was applied to reflect the degree of difference or correlation.

Results

Non-Canonical NF-κB Signaling Pathway Was Upregulated in LGG and Associated with Poor Patient Prognosis

In this study, we first examined the expression level of the core members (NF-κB1, NF-κB2, RELA, RELB and REL) of NF-κB signaling pathway. As shown in Figure 1A, these core members were universally upregulated in LGG relative to healthy tissues. Given that these core genes overlappingly involve in three NF-κB signaling branches: canonical NF-κB signaling, non-canonical NF-κB signaling and atypical NF-κB signaling, we further evaluated the activity in each branch using ssGSEA. Sample-wise correlation in Figure 1B showed that the expression of core genes highly correlated with activity of non-canonical NF-κB signaling (R=0.82) while only moderately correlated with canonical NF-κB signaling and atypical NF-κB signaling (R=0.57 and R=0.61, Supplementary Figure 1A and B), implying that the activation of NF-κB signaling is mainly attributed by its non-canonical branch. We also explored the association of NF-κB signaling with LGG patients’ survival. 410 LGG patients were stratified into low-, medium- and high non-canonical NF-κB signaling groups using X-Tile. Kaplan-Meier estimator showed that, among all three NF-κB signaling branches, only non-canonical NF-κB signaling stratified LGG patients into subgroups with differentiated survival outcome. (Figure 1C, p-value=0.016 and Supplementary Figure 1C and D). Specifically, LGG patients with high non-canonical NF-κB signaling activity own median survival time of 33.96 month, significantly shorter than LGG patients with medium or low NF-κB signaling activity. These results suggested that, compared with other branches, non-canonical NF-κB signaling imposes greater impact on the progression of LGG.

Figure 1 Activity of NF-κB signaling braches in LGG and their impact on survival outcome. (A) Box plots of normalized expression of core member genes of NF-κB signaling in LGG (red) and in healthy tissue (grey). Fold change and statistical significance were asterisked above (Mann–Whitney U-test, ***p<0.001). (B) Correlation between expression of core members of NF-κB signaling and activity of non-canonical NF-κB. (C) Kaplan-Meier curve of non-canonical NF-κB signaling.

Non-Canonical NF-κB Signaling Stratifies LGG Patients into Molecularly Distinct Groups

We thought to further characterize the molecular difference between LGG groups. Analysis on Whole Exome Sequencing (WES) reflected mutational spectrum in Figure 2A. We only displayed driver genes in LGG defined in Francisco Martínez-Jiménez et al’s study.²⁵ Non-canonical NF-κB signaling stratified groups were labeled at the top alone with major clinical information and genomic events in LGG. As expected, IDH1, ATRX and TP53 were prevalently mutated across all groups. EGFR was mutated in increasing trend from non-canonical NF-κB-low to -high groups (Fisher’s exact test, p-value < 0.01) and appeared in a mutually exclusive manner with IDH1 mutations. IDH mutations has been known as a prevalent genomic event correlated with favorable outcome of LGG patients.²⁶ On the other hand, mutations in EGFR may upregulate non-canonical NF-κB signaling and, along with other EGF induced oncogenic signaling pathways, promote tumor cell proliferation.²⁷ Thus, it is possible that the high frequency of EGFR mutants while low frequency of IDH1 mutation attribute to the poor survival of LGG patients in non-canonical NF-κB-high group. Besides, we found that CIC mutations and 1p19q co-deletion were enriched in NF-κB-low group (Fisher’s exact test, p-value < 0.01). Consistently, these two markers of oligodendroglioma highlighted the dominancy of oligodendroglioma in NF-κB-low group.²⁸

Figure 2 Genomic characteristics of LGG across non-canonical NF-κB groups. Mutation types were color marked. Major clinical information, genomic, epigenetic landmarks and pathological type were labeled on the top. Population frequency of each mutated gene was shown at left side.

Non-Canonical NF-κB Signaling Associated with Immunosuppressive Environment and Inhibition of Cytotoxic T Cells

TME mediates tumor progression, relapse and response to therapy^29,30 and NF-κB signaling has been known to modulate the TME to a great extent.^31,32 Therefore, we explored the association between activity of non-canonical NF-κB signaling and the landscape of Tumor-infiltrating immune cells (TIICs) in LGG samples. We first estimated the relative fraction of immune cell infiltration of each LGG sample. As Figure 3A shown, the activity of non-canonical NF-κB signaling positively correlated with overall degree of TIICs (Kruskal–Wallis test, p-value <0.05). This observation is consistent with our current knowledge of non-canonical NF-κB signaling in enhancing the immune response. TIICs represent a mixture of numerous immune cell types corresponding to different, even opposite, biological function. To gain a higher resolution of cellular composition of LGG samples under the influence of non-canonical NF-κB signaling, we further dissected TIIC subpopulations through deconvolution of RNA expression (See Method). We found that cytotoxic T cells demonstrated elevated infiltration from non-canonical NF-κB-low to -high group (Figure 3B). However, subsequent analysis found that the fraction of exhausted T cell showed the same trend (Figure 3C), implying that most cytotoxic T cells in NF-κB-high group may have lost their cytotoxic capacity against tumor cells. T cells exhaustion has been believed to be induced by some tumor cells capable of shaping its surrounding into immunosuppressive environment. Another frequently observed immune cell type with such capacity is macrophage type 2 (M2).We found that M2 was dominant across all three groups and showed elevated trend from non-canonical NF-κB-low to -high group (Figure 3D). Moreover, we analyzed correlation between the non-canonical NF-κB signaling and M2 markers. As shown in Figure 3E, the expression of M2 marker genes consistently correlated with activity of non-canonical NF-κB signaling. The above analysis results regarding cytotoxic T cells, exhausted T cells and macrophage were later validated in CGGA LGG dataset (Supplementary Figure 2). Together, we suggested that although the activity of nc NF-κB signaling positively correlates with immune infiltration, the elevated immune infiltration in nc NF-κB group displayed a detrimental immune microenvironment which may be a factor contributing to the poor prognosis of this group.

Figure 3 Immune cell infiltration of LGG samples. (A) Comparison of immune cell infiltration between non-canonical NF-κB signaling groups. (B–D) Comparison of cytotoxic T cells, exhausted T cells and macrophage infiltration between non-canonical NF-κB groups. (E and F) Heatmap of macrophage M2 marker gene expression and immune checkpoint expression in each non-canonical NF-κB signaling groups. (G–I) Comparison of TMB, IPS and GEP score between non-canonical NF-κB signaling groups. Statistical significance between groups was asterisked to the top (Kruskal–Wallis test; *p<0.05, **p<0.01 and ***p<0.001).

We hypothesized that the activation of the non-canonical NF-κB signaling pathway might be involved in inhibiting or exhausting cytotoxic T cells. The heatmap of the relative expression levels of T cell checkpoint, including PD-L1, PD-L2, TIM3, CTLA4 and SERPINB9 were elevated in the non-canonical NF-κB-high group (Figure 3F). The TIM3 and PD-L1 pathways have been identified as markers of T cell exhaustion.^33,34 These observations, alone with higher mutation burden in non-canonical NF-κB-high group (Figure 3G), implied that LGG patients of this group may be more likely to benefited from immunotherapy. We further assessed the immunogenicity in a more comprehensive manner. Immunophenoscore (IPS) reflects tumor’s potential responsiveness to immune checkpoint inhibitors (ICIs) by simultaneously measuring the expression of key genes involved in immune-effective cells, immunosuppressive cells, MHC complex and immunomodulators.²² As Figure 3H shown, LGG samples in non-canonical NF-κB-high group demonstrated significantly higher immunogenicity. Similarly, T cell–inflamed GEP is a rigorously validated model which measures the potential of antitumor response.³⁵ As expected, samples in non-canonical NF-κB-high group showed much higher GEP score among three groups (Figure 3I). In the past, clinical trials of immunotherapy in brain cancer have not been as successful as in non-CNS cancer partially due to anatomical structure and specialized immune environment of brain.³⁶ However, accurate molecular stratification may help identify LGG patients who are more likely to benefit from ICIs.

Non-Canonical NF-κB Signaling Pathway is an Independent Prognostic Factor for LGG Patients

Our previous analysis indicated that non-canonical NF-κB signaling significantly associated with worse OS of LGG patients (Figure 1C, Log rank test, p-value=0.016), we thought to explore the feasibility of constructing a simplified prognostic model that predict the likelihood of patient being alive after specific period of time (See Method). We randomly divided 410 patients in the TCGA-LGG dataset into the training set (n =273) and the testing set (n =137). Using training set, we prioritized nine feature genes for model construction through a series of feature selection technique. Figure 4A showed the ranked sample risk alone with the expression of prioritized feature genes. LGG patients were divided into low- and high-risk group based on the median risk score. The Kaplan-Meier curve showed significant survival difference with p-value less than 0.01 (Figure 4B and C). This survival difference was validated in testing dataset (Figure 4D and E). AUROC of the model was 0.85, 0.83 and 0.71 for 1, 3, and 5 years, respectively (Figure 4F). To test the independence of the model, potential covariates (tumor grade, age, gender, MGMT promoter methylation status, 1p19q status and IDH) were incorporated into Cox proportional-hazards model. Results indicated that, in addition to known survival predictors IDH and tumor grade, the risk model can independently predict survival of LGG patients (Figure 4G). Lastly, the prognostic independence of non-canonical NF-κB signaling as well as the robustness of the model was validated in an external validation set was retrieved from CGGA (Supplementary Figure 3).

Figure 4 Prognostic model of LGG. (A) Prognostic analysis of non-canonical NF-κB groups. LGG patients were divided into low- and high-risk groups based on median risk score. (B) Kaplan-Meier curve showed significant survival difference between low- and high risk groups. (C) ROC curve demonstrated the AUC for predicting patient status. (D–F) Prognostic performance during internal model validation using 137 LGG samples from TCGA. (G) Multivariate Cox model showed Independence of prognostic ability in predicting survival of LGG patients after adjusting known genomic and clinical confounders (*p<0.05, **p<0.01 and ***p<0.001).

Discussion

In this study, we identified that the non-canonical NF-κB signaling is upregulated in LGG and its activity is positively correlated with increased M2 macrophage infiltration and T cell exhaustion. More specifically, LGG patients with high non-canonical NF-κB activity display immune suppressive microenvironment, more complexed mutational landscape and poor survival outcome. We showed that activity of non-canonical NF-κB signaling is a validated prognostic signature in predicting the survival of LGG patients.

Our study showed that the non-canonical NF-κB signaling was upregulated in LGG and associated with poor OS rate of LGG patients. This result aligned with previous studies that non-canonical NF-κB promotes glioma cell invasion independently of the canonical NF-κB signaling.^37,38 Meanwhile, LGG patients in non-canonical NF-κB-high group displayed significantly higher frequency of EGFR mutation. Previous studies have demonstrated that constitutive activation of EGFR can trigger NF-κB signaling through activation of IκB kinase (IKK).²⁷ The EGF- dependent NF-κB activation, alone with other EGF induced oncogenic signaling, facilitates cancer cell proliferation and consequently resulted in poor survival of LGG patients. While patients in NF-κB-low and NF-κB-medium did not show survival difference, they differed in histology. Compared with NF-κB-medium group, NF-κB-low group was dominated by oligodendroglioma featured by high frequency of 1p19q co-deletion as well as CIC truncating mutations.

LGG is thought to be an immune desert relative to other cancer types. However, our analysis demonstrated that immune cells were infiltrated in LGG. In particular, the dominant infiltration of M2 macrophage was consistent with previous studies.^39,40 The fraction of M2 macrophages was the highest in the non-canonical NF-κB-high group. Additionally, the expression of M2 macrophage markers were remarkably upregulated in the non-canonical NF-κB-high group, suggesting that the activation of the non-canonical NF-κB signaling correlated with the immunosuppressive microenvironment in LGG.

T cells are critical in antitumor immunity. Non-canonical NF-κB signaling is essential for the generation and maintenance of effector T cells and memory T cells.¹⁶ We found that the activation of the non-canonical NF-κB signaling positively correlated with the overexpression of T cell checkpoint PD-L1, PD-L2, CTLA4, TIM3 and SERPIN9B. PD-L1, PD-L2 and CTLA4 are well-studied immune checkpoints that inhibit cytotoxic function of T cells. TIM3 functions as an immune checkpoint by prompting T cell exhaustion in human and mouse tumors.⁴¹ Previous studies demonstrated that cotargeting TIM3 and PD1 can induce tumor regression by reversing T cell exhaustion and restoring antitumor immunity in mouse colon carcinoma,⁴² Furthermore, researches have shown that TAMs are able to upregulate PD-L1 expression⁴² and targeting TAMs is one of the strategy under clinical investigation.^43,44 In our study, the molecular characterization suggested possible efficacy of those therapeutic approaches in treating LGG patient in non-canonical NF-κB-high group. Of course, further in-vitro experimental investigation need to be carried out to validate the therapeutic rationale.

Conclusion

Our results demonstrated that a prognostic signature based on the non-canonical NF-κB signaling is an independent prognostic factor for LGG patients. Two characteristics revealed by our study may explain the poor prognosis of LGG patients with elevated non-canonical NF-κB signaling: the distinct mutational landscape of this group may confer to the overcapacity of tumor cell proliferation; The immunosuppressive microenvironment featured with elevated M2 macrophage infiltration and cytotoxic T cell exhaustion may protect the tumor from immune response and further facilitate tumor growth. These molecular and cellular characteristics suggest potential response to immunotherapy of patients in this group. Meanwhile, our analysis demonstrated that, through feature selection technique, the prognostic model based on non-canonical NF-κB signaling can predict survival of LGG patients in a cost-effective way.