The immune regulation of BCL3 in glioblastoma with mutated IDH1

Background: Glioblastoma in the brain is the most malignant solid tumor with a poor prognosis. Screening critical targets and exploring underlying mechanisms will be a benefit for diagnoses and treatment. IDH1 mutation (R132) was used to distinguish glioblastoma grade and predict prognosis as a significant marker. However, the manner of IDH1 mutation regulating glioblastoma development was still unclear. Methods: To study the function of IDH1 mutation, multi-type sequencing data (transcriptome, methylation and copy number variation) from the GEO and TCGA database were analyzed using bioinformatics techniques. The biological functions of IDH1 mutation (R132) would be comprehensively evaluated from the regulatory networks, tumor immune microenvironment and clinical relevance. Then the analysis result would be validated by experimental techniques. Results: Compared with adjacent tissues, IDH1 was up-regulated in glioblastoma, which also positively correlated with the malignant degree and a poor prognosis. To further study the mechanism of mutated IDH1 (R132) function, 5 correlated genes (FABP5, C1RL, MIR155HG, CSTA and BCL3) were identified by different expression gene screening, enrichment analysis and network construction successively. Among them, the BCL3 was a transcription factor that may induce IDH1expression. Through calculating the correlation coefficient, it was found that in IDH1mut glioblastoma, the dendritic cell infiltration was reduced which may result in a better prognosis. In addition, the level of IDH1, FABP5, C1RL, MIR155HG, CSTA and BCL3 might also influence lymphocytes infiltration (eg. CD4+ T cell) and chemokine expression (CXCL family). Conclusions: IDH1 may participate in pathological mechanisms of glioblastoma via expression alteration or gene mutation. Furthermore, IDH1 mutation might improve prognosis via suppressing the expression of FABP5, C1RL, MIR155HG, CSTA and BCL3. Meanwhile, it was identified that BCL3 might perform similar immunomodulatory functions with IDH1 as an upstream transcript factor.


INTRODUCTION
Glioblastoma (GBM) is the most malignant and frequent cancer type in the central nervous system (CNS). Most GBM (~90%) develop rapidly in elderly patients. The unclear pathogenesis leads to therapeutic difficulty and high mortality. Despite the similar histologic appearance, the primary and secondary GBM are distinct cancers that originate from various precursor cells and require different therapeutic approaches. In 2008, the high-throughput sequencing of WHO grade I-IV gliomas were performed. In 12% of samples, a novel mutation at codon 132 (R132) of isocitrate dehydrogenase 1 (IDH1) was identified in ~80% of secondary GBM as a decisive genetic signpost [1][2][3][4]. Meanwhile, in primary GBM, this mutation was occurred rarely (<5%) [5,6]. Previous studies have identified that in secondary GBM, the IDH1 (R132) mutation is a more objective and reliable diagnostic marker than clinical observation [7,8]. The IDH1 (R132) mutation had never been linked to GBM before but now the mechanism of GBM development is under intense investigation.
B-cell CLL/lymphoma 3 (BCL-3), which belongs to the IκB family, can bind NF-κB homodimeric complexes of p50 or p52, which switches the transcriptional properties. BCL-3 overexpressed in breast cancer [21], nasopharyngeal carcinoma, endometrial cancer [22], hepatocellular carcinoma [23] and colorectal cancer [24] have been identified. Functionally, BCL-3 could participate in regulating the colony formation and cell cycle progression by regulating ubiquitination-mediated degradation of c-Myc in colorectal cancer [25]. Tu et al. [23] have reported that BCL3 promotes hepatocellular carcinoma growth by regulating cell proliferation and cell cycle through cyclin D1. However, the status and interaction with IDH1 of BCL3 in GBM have not been investigated.
In the present study, via bioinformatic analysis of GBM transcriptome sequencing datasets from GEO and TCGA databases and in vivo verification, try to find out the potential mechanism of IDH1 mutation of GBM development. Firstly, all samples were divided into wild type and IDH1-mutant groups. The different expression genes were identified according to |log FC|≥1.0 and p-value≤0.05 followed by enrichment analysis. In our study, via constructing a molecular network, the underlying mechanism of IDH1 mut -BCL3 in GBM development was explored. The results laid a foundation for further studying the function of IDH1 mutation, and also provided theoretical support for better screening markers/targets and designing the therapeutic schedule.

Sequencing datasets
On the one hand, the public GBM sequencing dataset (GSE122586) was obtained from the GEO database. According to the corresponding study [26], the GBM samples were collected in Beijing Tiantan Hospital from January 2005 to December 2009. The detailed sample information was introduced in Table 1. All surgically resected samples were rinsed with normal saline and then stored in liquid nitrogen immediately until use. All patients have voluntarily signed informed consent forms. The use of all human samples and the experimental procedures in that study were reviewed and approved by the Ethics Committee of Chinese Academy of Medical Sciences. In addition, another dataset GSE80729 contained U87 cells transfected with siRNA targeting BCL3 or control siRNA was used. The U87 glioma cell line is derived from a female patient with pleomorphic glioma. It differs from other glioma cell lines like U251 in cell proliferation, migration, and invasion. For example, U87 cells have exhibited a greater capacity for migration and invasion.
On the other hand, the RNA-seq data (RSEMnormalized), methylation array data (Illumina Human Methylation 450) and CNV data (Affymetrix SNP 6.0) of GBMs from the TCGA database were downloaded from the NIH National Cancer Institute GDC Data Portal (https://portal.gdc.cancer.gov/). Overall survival (OS) was identified from the diagnosis date until death or the end of follow-up. In addition, disease-free survival (DFS) was defined as the period from diagnosis until the first disease progression with the clinical sign.

Microarray data and enrichment analysis
Total RNA from cells and tissues was isolated using Trizol extractions (Invitrogen). The RNA quantity was Followed by background deletion, quantile normalization, and probe assembly. Different expression genes (DEGs) between normal vs. tumor tissues were detected by the empirical Bayes method [27]. The p-values were adjusted for multiple comparisons using the Benjamini-Hochberg procedure [28]. Genes with adjusted p-value < 0.05 and |log FC| ≥ 1.0 were considered as differentially expressed. Enrichment analysis of DEGs was performed with DAVID [29] and ClueGO [30]. The enriched GO (BP: biological process; CC: cellular component; MF: molecular function) and pathway terms were listed with participant genes [31]. Some other databases used are listed in Table 2.

Survival curves
Overall survival analyses were performed using the R package survival [32], and the patients were dichotomized based on the median expression. Kaplan-Meier estimator of survival was used to construct the survival curves. Log-rank tests (corresponding to a twosided z test) were used to compare overall survival between patients in different groups. The hazard ratio (HR) (95% confidence interval) was provided for comparison of the two groups. The p-values were adjusted for multiple testing based on the false discovery rate (FDR) according to the Benjamini-Hochberg method [33]. Proportional hazard assumptions were tested.

RT-qPCR
Total RNA was extracted using TRIzol® reagent (Takara Bio, Inc.). The concentration and purity of RNA were detected by NanoDrop®ND-1000 spectrophotometer. Total RNA was reverse transcribed into cDNA (10 µg) using the PrimeScript™ RT reagent kit (Takara Bio, Inc.), according to the manufacturer's protocol. The reaction conditions of the reverse transcription were: incubation at 25° C for 10 min, additional incubation at 42° C for 30 min, and heating at 95° C for 5 min. Then cDNA was diluted with DEPC water. Fluorescent RT-qPCR was performed according to the manufacturer's protocol. (Thermo Fisher Scientific, Inc.). The primers were designed and synthesized by Chongqing Life Biological Technology, Ltd. The reaction system volume was 25 µL consisting with cDNA (1 µL), 10X PCR buffer (2.5 µL), 10 mmol/l dNTPs (2 µL), PCR upstream primers (1 µL), PCR downstream primers (1 µL), Taq DNA polymerase (1 µL), and deionized water (16.5 µL). The reaction conditions were as follows: pre-denaturation at 95° C for 5 min, denaturation at 94° C for 1 min, annealing at 54° C for 45 sec, extension at 72° C for 1 min, in a total of 30 cycles, followed by extension at 72° C for 10 min. β-actin was used as the internal reference. The experiment was repeated 3 times independently.

Immunohistochemistry (IHC)
Before IHC staining, glioma specimens and normal brain tissues were fixed with 10% formalin and embedded with paraffin. The clinical samples were collected from clinical operation (Table 3). IHC staining followed standard protocol to evaluate the expression level of IDH1 and BCL3 in human tissues. Staining was graded on a scale of 0-3 according to the intensity and the percentage of immune-positive cells as follows: 0: no staining or <10% positive cells; 1: weak staining in >10% of cells or moderate staining in 10-70% of cells; 2: moderate staining in >70% of cells or strong staining in 10-70% of cells; 3: strong staining in >70% of cells.

Subcutaneous xenograft studies
All experimental mice were purchased from the Chongqing Medical University. To further study the role of IDH1 in GBM growth, 1 × 10 7 U87 cells with different treatments were injected into the axilla of 6 weeks old BALB/c nude mice. Tumor growth was determined by length (L) and width (W), which were measured at 4 weeks after the injection. The tumor volume (V), was calculated by V (mm 3 ) = (L×W 2 ) × 0.5. The use of all mouse samples and the experimental  procedures in this study were reviewed and approved by the Ethics Committee of Chongqing University Three Gorges Hospital.

Statistical analyses
Available samples from TCGA data were adequate because sufficient power using equivalent tests was observed in a previous study [34]. To test for differential expression across two groups (tumor and normal), the R package DESeq2 was used on raw count data [35]. The p-values were adjusted for multiple testing based on the false discovery rate (FDR) according to the Benjamini-Hochberg approach [36]. For comparison of two patient groups, the two-sided Student's t-test and Wilcoxon-rank sum test was used. For comparisons among multiple-patient groups, oneway ANOVA and Tukey's honest significant difference (HSD) post hoc tests were used. Distributions of data are shown either as individual data points, as box-andwhisker plots, or as violin plots. The p-values were adjusted according to the Benjamini-Hochberg method.

The expression of IDH1 in GBM
Abundant evidence suggested that IDH1 was functional in GBM progress while being regarded as a significant marker for the tumor classification. According to the RNA-seq of GBM from TCGA, it was detected that the IDH1 was up-regulated in GBM when compared with para-carcinoma tissue ( Figure 1A). In addition, the expression level might be positively correlated with malignant degree ( Figure 1B). The results may suggest that IDH1 plays a crucial role in tumorigenesis and malignant progression. The gene expression would be influenced by molecular characteristics like methylation, copy number and mutation. However, in GBM, the IDH1 expression was not associated with them. (Figure 1C-1E). Though the mechanism of IDH1 dysregulation was still unclear. The clinical data statistics showed that the lower IDH1 level could improve prognosis (overall survival not disease-free survival) ( Figure 1F, 1G). To further confirm the carcinogenesis of IDH1 in GBM, in the subcutaneous xenograft mouse model, once IDH1 expression was increased, the tumor growth rate was accelerated ( Figure 1H). Thus, it was speculated that once the IDH1 mutation appeared, the IDH1 enzyme activity would be lost and lead to GBM growth suppression and prognosis improvement.

The IDH1 mutation in GBM patients
Gene mutation often appears in cancer patients. Through analyzing whole-genome sequencing data of GBM patients from TCGA, the molecular aberrations of IDH1 were detected in nearly 7% of samples (Figure 2A). AGING  Among them, about half of patients appear the classic base mutation (R132H/G) ( Figure 2B). Statistical survival time showed that once IDH1 mutated, the prognosis will be better without significant difference (p-value=0.258) ( Figure 2C). In the previous section, it was demonstrated the expression and mutation of IDH1 were uncorrelated. To further confirm this conclusion, one transcriptome sequencing data (GSE122586) was detected and there was no significance between wild and mutated types ( Figure 2D). All results may suggest that the IDH1mutation may suppress GBM development by losing key enzyme activity, like isocitrate binding.

An enrichment analysis of DEGs regulated by IDH1mutation
In addition to altering their expression and functions, some mutations may also cause similar changes in downstream or correlated genes. To study the underlying regulatory mechanism of IDH1 mutation, the DEGs in IDH1 mut GBM were screened with |log FC|≥1.0 and p-value≤0.05 in both GEO and TCGA datasets ( Figure 3A, 3B). To ensure accuracy and increase credibility, the dysregulated gene screened in both databases (n=102) were used for the next analysis ( Figure 3C, 3D). Furthermore, there are 6168 genes associated with IDH1 in GBM. Among them, only 16 genes (down-regulated: PDPN, TUBA1C, LGALS3,  FABP5, ANXA2P1, ANXA2, C1RL, KCNE4, LGALS1, DCDC2, MIR155HG, CSTA, BCL3 and FBXO17; up-regulated: RANBP17 and DLL1) were dysregulated in IDH1 mut GBM simultaneously ( Figure 3E and Supplementary Figure 1). Following enrichment analysis for GO function (BP, MF and CC), some immune functions like the immune effector process (adjust p-value=7.831×10 -6 ) and immune system development (adjust p-value=1.711×10 -4 ) were enriched. It was hypothesized that IDH1 mutations could mediate tumor development via regulating the immune microenvironment.

The immune correlation of IDH1 mutation in GBM
In most cancer studies, the immune microenvironment was often considered as an extracellular factor affecting tumor development. Tumor cells can achieve immune escape or suppression by dysregulation of specific genes. It has been demonstrated that IDH1 mutation could be involved in immune functions by regulating the expression of multiple genes. To study the impact of IDH1 mutation on the GBM immune microenvironment more comprehensively, the correlation between IDH1 mutation and lymphocyte infiltration, immune inhibitor, immunostimulator, MHC molecule, chemokine and chemokine receptor was calculated ( Figure 4A). Almost lymphocyte infiltrations (Th1 cell) and all MHC molecule (HLA family) expression were negatively correlated with IDH1 mutation. Meanwhile, the expression of some immunostimulators (eg. KLRK1) and chemokine (eg. CCL3) were positively correlated with IDH1 mutation. Lymphocyte infiltration is the most critical event directly regulating tumor immune escape or suppression in the microenvironment. In IDH1 mutated samples, the infiltration of B cell (p-value=0.018) and CD4+ T cell (p-value=0.016) were promoted while CD8+ T cell (p-value=0.023) and dendritic cell (p-value=0.017) were suppressed ( Figure 4B). In addition, it was verified that only the infiltration of dendritic cells could significantly alter prognosis (Log-rank P=0.002).

The clinical character of DEGs correlated with IDH1
As mentioned above, a total of 16 DEGs were identified in IDH1mutated samples which also correlated with IDH1 expression. Among them, only 5 positively related genes (FABP5, C1RL, MIR155HG, CSTA and BCL3) would significantly induce the poor prognosis in GBM patients (p-value≤0.05) ( Figure 5A). Notably, in the enrichment analysis above, FABP5, C1RL, CSTA and BCL3 could also participate in various immunology processes, like mature B cell differentiation. Therefore, it was concluded that FABP5, C1RL, MIR155HG, CSTA and BCL3 might mediate the functions of IDH1 mutation in the GBM immune microenvironment. The correlations between the gene expression and lymphocyte infiltrations in the GBM patients were calculated (Supplementary Figure 2). Among them, MIR155HG, a non-coding gene, was not correlated with all types of lymphocyte infiltrations (p-value≥0.05). However, the dendritic cell infiltration was most significantly positively correlated with other genes level (p-value≤0.05). Compared with normal tissues, FABP5, C1RL, MIR155HG, CSTA and BCL3 were significantly increased in GBM ( Figure 5B). Nevertheless, in IDH1 mut GBM samples, all 5 genes were down-regulated significantly ( Figure 5C). Besides MIR155HG, other genes also own higher molecular aberrations frequency in GBM ( Figure 5D). It may suggest that in GBM, FABP5, C1RL, CSTA and BCL3 could be regarded as oncogenes like IDH1.

BCL3 may regulate IDH1 as a transcription factor
The genes screened previously, including IDH1, perform functions in an interactive network rather than independent manners. Based on the STRING database, an interactive network was constructed for systematically understanding gene functions ( Figure  6A). In this network, BCL3 may interact with most genes, including SP1, TP53 and MYC which were defined as crucial oncogenes. Moreover, IDH1 was predicted as a target of BCL3. In GBM tissue, BCL3 was increased and positively correlated with IDH1 ( Figure 6B and Supplementary Figure 2). To further confirm that BCL3 could target IDH1, possible action sites were screened by ChIP-sequencing datasets. A total of 4 predicted sites were verified on IDH1 gene, among which one was in the promoter region near the TSS (transcription start site) (Chr2, 208265841-208266367, signal value=4.69), others were in the gene body. The maximum signal peak was identified at Chr2, 208255168-208255438 (signal value=5.50) where has a higher GC percent ( Figure 6C). However, whether BCL3 regulates IDH1 and the specific binding site needs to be further studied.
To confirm that BCL3 can positively regulate IDH1 in glioma, transcriptome sequencing data containing the glioma cells transfected with siRNA targeting BCL3 or control siRNA were analyzed (Figure 7). After DEG screening, a total of 355 up-regulated genes and 140 down-regulated genes. In addition, the expression of IDH1was downregulated when BCL3 expression was disrupted. The result may further be confirmed that the BCL3 is a transcription factor of IDH1.

BCL3 and IDH1 in GBM immune microenvironment
IDH1 mutation would alter the expression level of BCL3 ( Figures 5C, 6D). Meanwhile, BCL3 was predicted to regulate IDH1 as a transcription factor. In addition, both BCL3 and IDH1 might participate in immune functions, like dendritic cell infiltration. To comprehensively study the mechanism of IDH1 and BCL3 in tumor immunity regulation, the correlation between the expression level and various immune events was calculated ( Figure 8A). Interestingly, the result showed that BCL3 and IDH1 own a similar correlation of lymphocyte infiltration and immune factor expression. The correlated degrees of BCL3 were higher than those of IDH1 maybe also indicate that in the immune microenvironment, BCL3 owns a larger regulatory capacity. In different subtypes, the similar expression tendency of BCL3 and IDH1 further confirms they have a close relationship.

DISCUSSION
In a population-based study, the median overall survival of clinically diagnosed GBM patients with IDH1 mut was 234 days, significantly longer than IDH1 wt patients (141 days; p-value=0.003) [37,38]. Furthermore, an analysis of GBM patients treated with surgery or/and radiotherapy showed that the mean overall survival of IDH1 mut patients was 813 days, longer than patients with IDH1 wt GBM (339 days; p-value <0.0001) [39]. However, the underlying mechanism of IDH1 mutation regulating the prognosis of GBM patients was still indistinct. In the present study, via bioinformatics analysis of multiple level sequencing datasets from GEO and TCGA databases, it was demonstrated that IDH1was up-regulated in GBM when compared with adjacent tissue. And the IDH1 mutation did not significantly alter its expression level. However, in Bleeker's study [15], it was found that NADPH production is hampered in GBM with IDH1 (R132) mutation. Moreover, mutated IDH1 consumes rather than produces NADPH. Thus, likely reducing NADPH levels even further. The low NADPH levels may sensitize GBM to irradiation and chemotherapy, thus explaining the prolonged survival of patients with mutated glioblastoma [40,41].
Through analysis of transcriptome sequencing data, FABP5, C1RL, MIR155HG, CSTA and BCL3 were screened which were positively correlated with IDH1 and down-regulated in IDH1 mut GBMs. However, whether mutated IDH1 regulate these genes expression and the mechanism is still unknown. According to the references [15], once IDH1 mutated, the function would be changed. The tumor metabolite 2-HG could be produced from α-KG by accelerating NADPH consumption [13,14]. The increased 2-HG competitively suppresses the α-KGdependent dioxygenases which mediate DNA demethylation [16]. Thus, in GBM, it was speculated that IDH mutation regulates gene expression in an epigenetic manner like DNA hypermethylation rather than direct regulation.
A molecular network including IDH1 and correlated genes was constructed. Among them, lower levels of FABP5, C1RL, MIR155HG, CSTA and BCL3 could improve the prognosis in GBM patients. To understand the function of this network, the enrichment analysis showed that the critical genes (FABP5, C1RL, CSTA and BCL3) could also participate in various immunology processes. It has been studied that in GBM, FABP5 expression correlated with an undifferentiated tumor phenotype as a known tumor-associated antigen that could respond to B cell [42,43]. It was also confirmed that C1RL was upregulated in GBM, especially mesenchymal GBM and primary GBM. Increased C1RL expression accompanied the IDH1-wt phenotype in both lower-grade glioma (LGG) and GBM. C1RL-associated genes were mainly enriched in biological processes related to the immune response. C1RL expression was also correlated with reduced tumor purity and increased M2 macrophage infiltration [44]. Followed by the immune score, it was hypothesized that the IDH1 may combine with BCL3 to change the GBM immune microenvironment like dendritic cell infiltration. In addition, CSTA upregulation has previously been described in human malignant gliomas: CSTA positive cells in GBM tumor samples were located close to tumor blood vessels, particularly in leukocytes and inflammatory host cells, possibly reflecting the level of inflammatory cells in the tumor tissue. CSTA expression displayed a significant correlation with markers (CD68 and CXCR4) of invasive GBMs [45]. BCL3 was previously found to be a putative protooncogene in human cancers and attenuates the efficacy of temozolomide in glioblastoma cells [46]. Only one published study has proved that in the mouse xenograft model, BCL3 inhibited tumor growth via positively regulating STAT3/p-STAT3 and the downstream targets including cyclin D1 [47]. In our study, according to ChIP-sequencing data, multiple BCL3 predicted sites on the IDH1 gene sequence were identified. Therefore, it was speculated that BCL3 might regulate the expression of IDH1 which also needs further research. Through functional enrichment analysis, both IDH1 and BCL3 were found to be involved in multiple immune activities. Previous studies have found that BCL3 could be involved in immune response in various cancer types [48][49][50].
After calculating immune correlation, lots of similarities were found between both 2 genes. For example, both IDH1 (partial. cor=0.289) and BCL3 (partial. cor=0.512) were positively correlated with dendritic cell infiltration. Dendritic cells are the most powerful antigen-presenting cells. DC-based vaccination has the potential to target and eliminate GBM cells and enhance the responses of malignant cells to the existing therapies with minimal damage to healthy tissues [51]. The efficacy of this therapy can be strengthened in several ways like regulation of regulators in the GBM microenvironment [52]. This will provide a basis for better research on the role of IDH1 mutation in GBM, and has guiding significance for the development of new therapeutic protocols based on the immune function of IDH1 and BCL3.

Availability of data and material
All data and material were available. The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

Ethics approval and consent to participate
The ethical approval about all animal experiment and clinical samples used in our study was obtained from the Chongqing University Three Gorges Hospital.

AUTHOR CONTRIBUTIONS
Shibing Fan, Na Wu, Shichuan Chang, Long Chen and Xiaochuan Sun conceived the project. Shichuan Chang, Na Wu and Shibing Fan analyzed the experiment data. The animal experiment was performed by Shibing Fan and Na Wu. Long Chen and Shichuan Chang acquired sequencing datasets and interpreted the results. Shibing Fan and Long Chen and wrote the manuscript. Xiaochuan Sun polished up this article.

ACKNOWLEDGMENTS
The results shown here are in part based upon data generated by TCGA Research Network (https://cancergenome.nih.gov). The authors would like to acknowledge the helpful comments on this paper received from the reviewers.