The High Expression of CD276/HAVCR2 and CD163 is an Adverse Immune Subtype of Glioblastoma and is Closely Related to Epithelial-Mesenchymal Transition

Background Glioblastoma (GBM) is one of the most malignant tumors that can aict the central nervous system. Previous studies have observed that there are individual differences in the treatment response of immune checkpoint inhibitors in glioblastoma. This study’s aim is to ascertain the factors that may affect the ecacy of immunosuppressant therapy. Methods The clinical data of this study were obtained from a public database. Then, the data was analyzed and processed by R software and corresponding R package. To verify the results of the analysis, information was gathered from 89 GBM patients in our hospital and thereafter the corresponding paran sections were stained and quantitatively analyzed by immunohistochemistry. Results From the analysis, it was observed that both CD276 and HAVCR2 were signicantly overexpressed in GBM and could be associated with patient prognosis. The analysis of single cell RNA sequencing data and GBM data analysis found an immune subtype with poor prognosis. Further analysis found that the high expression of CD276, HAVCR2 and CD163 was closely related to epithelial-mesenchymal transition (EMT) and could affect the patient prognosis of PD-L1 high expression. GSVA enrichment analysis showed that CD276, HAVCR2 and CD163 might induce EMT by JAK-STAT3 signaling pathway, and RUNX1 and IKZF1 might be transcription factors that regulate CD276/HAVCR2 high expression. Conclusions We found an immune subtype with poor prognosis of GBM, the high expression of CD276, HAVCR2 and CD163 with EMT are closely related and may be one of the factors affecting the ecacy of Anti-PD-L1. Family Zinc


Abstract
Background Glioblastoma (GBM) is one of the most malignant tumors that can a ict the central nervous system. Previous studies have observed that there are individual differences in the treatment response of immune checkpoint inhibitors in glioblastoma. This study's aim is to ascertain the factors that may affect the e cacy of immunosuppressant therapy. Methods The clinical data of this study were obtained from a public database. Then, the data was analyzed and processed by R software and corresponding R package. To verify the results of the analysis, information was gathered from 89 GBM patients in our hospital and thereafter the corresponding para n sections were stained and quantitatively analyzed by immunohistochemistry. Results From the analysis, it was observed that both CD276 and HAVCR2 were signi cantly overexpressed in GBM and could be associated with patient prognosis. The analysis of single cell RNA sequencing data and GBM data analysis found an immune subtype with poor prognosis.
Further analysis found that the high expression of CD276, HAVCR2 and CD163 was closely related to epithelial-mesenchymal transition (EMT) and could affect the patient prognosis of PD-L1 high expression. GSVA enrichment analysis showed that CD276, HAVCR2 and CD163 might induce EMT by JAK-STAT3 signaling pathway, and RUNX1 and IKZF1 might be transcription factors that regulate CD276/HAVCR2 high expression. Conclusions We found an immune subtype with poor prognosis of GBM, the high expression of CD276, HAVCR2 and CD163 with EMT are closely related and may be one of the factors affecting the e cacy of Anti-PD-L1. Background Glioblastoma (GBM) is a malignant tumor of the central nervous system. The main treatment methods are surgery, adjuvant radiotherapy and chemotherapy. Despite the active treatment, GBM exhibits characteristics such as short median survival, high recurrence rate, poor prognosis, etc. [1][2][3][4]. The primary method of adjuvant treatment after operation is to administer chemotherapy. However, the individual effect of chemotherapy is poor [5][6][7]. The low e cacy of traditional tumor treatment programs have prompted researchers to actively seek more effective and targeted treatment options. After active exploration, the researchers found a research direction that may change the current therapeutic effect: immune checkpoint inhibitors. However, the therapeutic effect is still not satisfactory to us, and it is critical to nd the potential factors that affect the curative effect.
Normally, when the body has an in ammatory response, the main function of the immune checkpoint is to reduce the molecules involved in the in ammatory response in the body and maintain immune homeostasis, thereby preventing the autoimmune response. But the tumor cells use this mechanism to achieve immune escape [8]. Therefore, the study of immune checkpoint is crucial. Immune checkpoint inhibitors primarily activate the body's immune response by inhibiting immune checkpoint gene surface ligands, thus affecting the development and progression of tumors. At present, the most studied immune checkpoint in GBM is PD-1 or PD-L1 [9][10][11][12], and the current GBM immune checkpoint inhibitors, mainly targeting PD-1 or its ligand, PD-L1 [12][13][14], have been used in clinic to improve the prognosis of GBM patients [15]. Some studies have shown that for GBM patients, the combination of immune checkpoint inhibitors and other drugs is more effective than treatment alone [16][17][18][19][20]. But previous studies have also observed that there are individual differences in the treatment of immune checkpoints and thus is effective only for some patients [11], and studies of other immune checkpoints have been poorly reported.
This indicates the need to dispense personalized treatment as per a patient's condition during the treatment process. Hence, it is necessary to explore the possible causes of individual differences in treatment and other effective immune checkpoint genes and subtypes with good prognosis, so that the personalized treatment of GBM is made possible.
Epithelial-mesenchymal transition (EMT) plays an important role in tumor invasion, metastasis and so on, and is closely related to its related inducers, transcription factors, signaling pathway genes and tumorassociated macrophages (TAMS), and TAMS can induce EMT [21,22]. Most of the immune cells in brain tumors are macrophages. TAMS are considered to be supportive matrix of tumor growth and are generally considered to be associated with CD163, CD206, etc. [23,24]. Most TAMS are considered to exhibit a M2 phenotype, whereas M2 macrophages can produce cytokines such as IL-10, IL-6 and TGF-β, which contributes to immunosuppression and thus facilitates tumor growth and induce EMT [22,24,25]. CD163 is primarily expressed by monocytes/macrophages and is a marker of M2-type macrophages. In tissues, CD163 is mainly resident in tissue macrophages [26], and CD163 is considered the most speci c marker of TAMS [27,28]. Recent studies have reported that CD163 can promote the occurrence of glioma through CK2 and may be a potential new target for glioma treatment [29,30], and closely related to EMT in other cancers [31,32]. But the association and function of CD163 with immune checkpoint genes and EMT in GBM is yet to be examined.
For the immune checkpoint genes in this study, previous researches only found CD276 to be associated with GBM prognosis [33,34]. HAVCR2 was also found to be associated with adaptive immune cells, which can affect in ammatory response and also affect the prognosis of GBM [17,35,36]. This study combined multi-database to study the relationship of them, thus ultimately identifying the high expression of CD276, HAVCR2 and CD163 may be one of the factors affecting the e cacy of Anti-PD-L1 by induce epithelial-mesenchymal transition.

Data source
In this study, the RNA expression data and clinical data of GBM patients was obtained from The Cancer Genome Atlas (TCGA) (n=168) (http://www.tcga.org/) database, Chinese Glioma Genome Atlas (CGGA) (n=138) (http://www.cgga.org.cn/) database, and single cell RNA sequencing data of GBM patients from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database: GSE84465(4 patient samples), GSE 103224(8 patient samples). While the immune and stromal scores were acquired from the ESTIMATE algorithm, the expression data of immune cells were obtained from TIMER database (https://cistrome.shinyapps.io/timer/). The pathway data for GSVA enrichment analysis was taken from the MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb/). To verify the analysis results of this study, 89 cases of pathological results of the past 5 years were collected from the a liated hospital of the Zunyi medical university. This consisted of high-grade glioma patients clinical information and corresponding para n tissue sections. All patients were diagnosed by two independent neuropathologists as per the World Health Organization classi cation guidelines. This study was approved by the ethics committee of the a liated hospital of the Zunyi medical university. To nd upstream transcription factors (TFs) that may regulate CD276/HAVCR2, we jointly UCSC (https://xenabrowser.net/datapages/), JASPAR (http://jaspar.genereg.net) and CISTROME (http://cistrome.org/) databases for analysis.

Data processing
Patients whose survival data was missing were excluded from this study. This reduced the dimensionality of the data of numerical variables to two categories of data based on the best cut-off value. Initially the heatmap analysis and the survival curve of immune checkpoint genes were undertaken. The expression of CD276 and HAVCR2 in each grade of glioma was analyzed in GEPIA (http://gepia.cancer-pku.cn/) database. To further determine the expression of genes in macrophages, TSNE clustering and data processing were performed on GSE84465 and GSE103224 dataset. Simultaneously, the survival curve and in ltration was plotted through the expression data of immune cells in the TIMER database to accurately comprehend the effect of immune cells on the prognosis of GBM patients. The correlation between CD276/HAVCR2 and CD163 was further assessed through the application of the Pearson correlation analysis. Then CD163, CD276 and HAVCR2 were regrouped as per the results through the cut-off best truncation values. To grasp the possible biological processes involved in TAMS, CD276 and HAVCR2, the GSVA enrichment analysis was performed.

EMT related genes analysis and lasso risk model construction
We analyzed the relationship between I-IV group and EMT related genes by T test and analysis of Variance. We further calculated the EMT risk score of each patient by using the Lasso algorithm using genes related to the prognosis of patients, and then divided the patients into high or low risk group.

Immunohistochemistry
Immunohistochemistry was utilized to detect the protein expression levels of CD163, CD276 and HAVCR2. We collected the corresponding para n pathological sections of the patients. Anti-CD163, CD276 and HAVCR2 antibodies were used at a dilution of 1:200. Each stained section was reviewed and evaluated by two independent pathologists. The expression of the protein was quanti ed by Image-pro Plus 6.0 software. The quantitative formula was protein expression = IOD/Area.

Statistical analysis
All statistical analyses were undertaken using the R software (3.6.1). Statistical methods such as Cox regression analysis, Univariate variance analysis, Pearson correlation analysis, Log-rank test, T test Analysis of Variance and Logistic regression analysis were applied. For all statistical methods, P <0.05 was considered statistically signi cant.

Identi cation of immune checkpoint genes CD276 and HAVCR2
To identify the immune checkpoint genes of this study's focus, it was observed that CD276, HAVCR2 and CD47 genes were signi cantly overexpressed in most tumors, including GBM (Fig. 1a). It was noticed that CD276, CD47 and HAVCR2 were also highly expressed in the GBM samples of TCGA and GEPIA database ( Fig. 1b-d). Additional survival analysis of immune checkpoint gene revealed that CD276, TNFSF14, VTCN1, IDO1, CD40, TNFSF18, TIGIT, CD274, CD70 and HAVCR2 could in uence the prognosis of GBM patients (Fig. 1e-h and Figure s1). In summary, this study identi ed CD276 and HAVCR2 as the genes of interest.

Correlation analysis between CD276, HAVCR2 and macrophages
Macrophages are important TAMS cells which play an important role in GBM, especially in the M2 macrophages. It was observed that primarily the CD276, HAVCR2 had higher expression in macrophages after clustering RNA single-cell sequencing data, whereas its expression in the neuronal cells was low GBM poor prognostic immune subtype determination and the relationship with IDH status and clinical information While a positive correlation was observed between CD276, HAVCR2 and CD163, no signi cant correlation was noticed with NOS2 ( Fig. 3). By further re-grouping the data by the cut-off best truncation values, namely I (CD276 low CD163 low /HAVCR2 low CD163 low ), II (CD276 low CD163 high /HAVCR2 low CD163 high ), III (CD276 high CD163 low /HAVCR2 high CD163 low ), IV (CD276 high CD163 high /HAVCR2 high CD163 high ), it was found that the IV group had the worst prognosis, I group had the best prognosis, and the results were statistically signi cant (Fig. 4). Further analysis found that most of the IDH status in the IV group was wild type, and there was no signi cant relationship between gender and age (Fig. 5), which might imply that patients with group IV have higher potential for immunotherapy and have a better curative effect [37]. To verify this study's results, IHC staining was performed on para n sections of patients who had been in our hospital for the past ve years and were pathologically diagnosed as high-grade glioma (Fig.  6a). It was noted that the degree of staining of para n sections was different in different patients undergoing the same treatment. The patients were also grouped as per the same best cut-off value after the quanti cation of the sections. The results showed that patients with IV group had the lowest survival whereas I group had the highest survival, which was also statistically signi cant (Fig. 6b and c). The four groups were associated with immune score and stromal score: I group had the lowest immune score and stromal score, while IV group had the highest immune score and stromal score (Figure s5), and macrophage in ltration in IV group was the highest (Figure s6).

GSVA enrichment analysis
Through the analysis of the GSE84465 data set, it was noted that the macrophages were divided into 0,2 clusters, wherein 0 cluster were basically tumor-associated macrophages and the 2 cluster were basically paracancerous cells. Furthermore, through the GSVA enrichment analysis and differential expression analysis, it was ascertained that there were differential signaling pathways (Fig. 7a). After further calculation gene set score, it was associated with EMT and JAK-STAT3 signaling pathway (Fig. 7b and   c). GSVA analysis of CD276, HAVCR2 and CD163 showed that was positively correlated with macrophage activation, differentiation, EMT and JAK-STAT signaling pathway. This result can be veri ed in both TCGA and CGGA databases (Fig. 8).
High expression of CD276, HAVCR2 and CD163 was closely related to EMT and affected the prognosis of patients with high expression of PD-L1.
TAMS and EMT were closely related and enrichment analysis found signi cant enrichment in EMT signaling pathways, GSVA enrichment analysis of CD276, HAVCR2 and CD163 also found signi cant correlation with the EMT. Further analysis found that EMT most of the associated genes were signi cantly high expression in IV group (Fig. 9). Interestingly, IV group accounted for more than 30% of the patients in the PD-L1 high expression group, and IV group had a worse prognosis than the other three groups (Fig. 10). In order to further explore the reasons for the difference in prognosis, we establishing a risk model by EMT related genes, calculated the possible risk values for each patient (Figure s7), and found that IV group was also closely related to the high expression of EMT related genes in the PD-L1 high expression group, and were more likely to develop EMT high risk scores, and IV_EMT_score high had a worse prognosis (Fig. 11). We speculated that this might be one of the reasons for the individualized differences in Anti-PDL1 treatment.

Determination of upstream transcription factors
We analyzed the differential expression of tumor-associated transcription factors (n=318) downloaded from the CISTROME database using TCGA GBM data combined with GET data, the screening conditions were |LogFC|>1, P<0.05, and nally 58 differentially expressed transcription factors were obtained. We further analyzed the correlation with CD276/HAVCR2, the screening conditions were |R|>0.5, P<0.05. We found that the positive correlation between RUNX1/IKZF1 with CD276/HAVCR2 was the highest (Fig.  12a), and the correlation coe cient was 0.7, 0.79, respectively. Meanwhile, we predicted the transcription factors that might bind to CD276/HAVCR2 by UCSC database (Minimum score=400), and then cross with CD276/HAVCR2 related transcription factors, nally, the RUNX1 and IKZF1 was obtained (Fig. 12b and c).
Further analysis found that The RUNX1/IKZF1 binding peaks were found in CD276/HAVCR2 promoter region by UCSC analysis (Fig. 12C and E). To search for speci c possible binding sites, we screened the possible binding sites in the CD276/HAVCR2 sense sequence or antisense strand by JASPAR database, the screening conditions were Relative pro le score threshold>90%, and nally the possible binding sites were tcctgcggttg, caaagaggaagc, respectively ( Fig. 12d and f).

Discussion
Immune checkpoint inhibitors are yet to attain a clear clinical bene ts in GBM [38]. Hence, an improved understanding of the possible causes of individual differences is needed to identify a subset of patients that may bene t from such treatments, thereby improving their current treatment. CD276 and HAVCR2 are key immune checkpoints in GBM, and its high expression is associated with poor prognosis in the patients [34,35,[39][40][41]. This observation is consistent with the results of this study. Furthermore, CD163 is considered the most speci c marker of TAMS, involved in the occurrence, development, migration, angiogenesis, EMT, etc., of GBM [42,43]. PD-L1 is the focus of current immunotherapy research, the study found that its expression is closely related to EMT [44,45], and EMT is the focus of current cancer resistance research [46,47]. But no research has been done to explore the causes of individualization of patient therapy through their relationship. So, an important purpose of this study to nd the additional factors that may affect the e cacy of immunosuppressant therapy and also consisted of a simple exploration of the mechanisms that may be involved.
Macrophages are important cellular components in tumor microenvironment. This study's cluster analysis of single-cell RNA sequencing data showed that CD276 and HAVCR2 was mainly highly expressed in macrophages and not in neuronal cells. Additionally, the TCGA and CGGA databases analysis showed that while CD276 and HAVCR2 had positive correlation with CD163, it had no signi cant correlation with NOS2. Results indicated that the existence of some close relationship between CD276, HAVCR2 and CD163. Hence, to further investigate the relationship between them, CD163, CD276 and HAVCR2 were re-grouped by the cut-off best truncation values. The outcomes showed that the patients with CD276 high/low expression also had CD163 high/low expression. The same phenomenon was exhibited in HAVCR2. The prognosis of CD163 high/low expression is different in patients with CD276/HAVCR2 high/low expression, and the prognosis of IV group is the worst, while that of I group is the best. This result was validated in our clinical data collection. Combined with IDH status analysis of patients, which might imply that patients with group IV have higher potential for immunotherapy and have a better curative effect [37]. The most widely studied immune checkpoint in the eld of cancer is PD-L1/PD-1. But study found that the therapeutic effect was not satisfactory and there were differences in therapeutic effect [48]. Current clinical studies of non-small cell lung cancer have assumed that people with PD-L1 high levels of expression bene t more from the use of immune checkpoint inhibitors than patients with low or no PD-L1 expression [49]. This prompts us to further explore in GBM according to this phenomenon, thus enabling us to further individualize treatment. We subsequent analysis found that CD276, HAVCR2, CD163 and EMT might closely related, so we analyzed the EMT related genes and found a signi cant positive correlation with the IV group. Interestingly the high expression group of PD-L1 included more than 30% IV group patients, the prognosis of the IV group was signi cantly worse than that of the other three groups. The IV group of PD-L1 high expression was also signi cantly positively correlated with EMT related genes and patients were more likely to develop EMT high risk scores. This suggests that the EMT molecular subtypes in PDL1 high expression group were mainly the IV group, and IV_EMT_score high had a worse prognosis. Therefore, the results showed that high expression of CD276, HAVCR2 and CD163 was closely related to EMT and affected the prognosis of patients with high expression of PD-L1, and might be one of the factors affecting the e cacy.
The possible involvement of them in biological processes and possible upstream mechanisms was also brie y explored in this study. The GSVA enrichment analysis and difference analysis of the GSE84465 data set highlighted the differences in several signaling pathways related to TAMS, including JAK-STAT, EMT and other signaling pathways. In addition, the CD276/HAVCR2 GSVA enrichment analysis also showed positive correlation with the positive regulation of EMT and JAK-STAT signaling pathway. This result can be veri ed in both TCGA and CGGA databases. Previous studies have shown a clear relationship between JAK-STAT signaling pathways, EMT and drug resistance [50], while TAMS can induce the occurrence of EMT [21], and studies in other tumors have found CD276,HAVCR2 and CD163 associated with AKT signaling pathways and EMT [51][52][53], but little is known about the relationship in GBM. To explore upstream regulatory factors that may affect CD276/HAVCR2 high expression in GBM. We analyzed tumor-associated transcription factors and found RUNX1/IKZF1 signi cantly positive correlation with CD276/HAVCR2, and corresponding binding peaks in the promoter regions of the CD276 sense sequence and the HAVCR2 antisense strand, indicating that there may be regulatory links between them. Previous studies have found that RUNX1/IKZF1 can induce the occurrence of EMT, thus affecting the occurrence and development of tumors [54][55]. But the relationship with CD276/HAVCR2 is not reported. we need further study to verify.
Nowadays, personalized treatment of patients is gaining in importance. Early detection of the key causes of individual differences in treatment of patients and the early formulation of personalized treatment can enable the patients to access to more treatment options. Though this study has reached the same conclusion through multi-database data analysis combined with the retrospective data that was collected in this study, the number of patients reviewed was small. Hence, more retrospective and prospective research support is required.

Conclusion
Despite the poor prognosis of GBM patients and the high postoperative recurrence rate and the use of immune checkpoint inhibitors provides additional treatment options for improving the prognosis of GBM patients. Through combined multi-database analysis, this study was able to identify immune subtypes with a good prognosis and factor that may affect the e cacy of Anti-PD-L1, thus providing a potential possibility for individualized treatment of immunity.       We show the sequencing analysis and quality control process RNA single cells by GSE84465 data sets.   (a-d) Scatter plots of CD276/HAVCR2 and CD163 in TCGA and CGGA databases, respectively, a linear regression line is plotted with the gray shaded region showing the 95% con dence interval. Pearson's correlation coe cient r and P values are given at the top. The median of CD276/HAVCR2 and cd163 expression is indicated by a gray dotted line. Log-rank statistics were applied to identify the optimal cutoff for transforming the continuous variable of gene expression into categorical high-and low-expression groups. The test score at each candidate cut-off across the log-transformed gene expression values was plotted. The highest value (indicated with a blue arrow) was using solid blue lines to best separating patients into 4 groups (named groups I to IV). (E-H) Represents the survival curves of I to IV groups in TCGA and CGGA databases, respectively. The prognosis was the best in group I and the worst in group IV (p 0.05).

Figure 5
In TCGA and CGGA database, the IDH status of IV group was mostly wild type, and there was no signi cant relationship between gender and age.     IV group accounted for more than 30% of the PD-L1 high expression group, and the prognosis of IV group was signi cantly worse than that of the other three groups, results are veri ed in the CGGA database.

Figure 11
We found that IV group of patients were more likely to have high EMT scores, it was also closely related to the high expression of EMT related genes in the PD-L1 high expression group, and IV_EMT_scorehigh had a worse prognosis, and the IV_EMT_scorelow prognosis was similar to that of I_III group, which also explains the reason of difference in survival between groups in PD-L1 high-expression group. * 0.05, * * 0.01, * * * 0.001, **** 0.0001.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.