Immune response related genes are found to be up-regulated in EGFR-MUT LGGs.
To explore the characteristics of EGFR mutations in lower grade glioma, we performed differential analysis of gene expression using RNA-seq data from TCGA-LGG cases with and without EGFR mutations. A total of 499 LGG cases, 471 of which had wild-type-EGFR (EGFR-WT) and 28 had mutant-EGFR (EGFR-MUT), were thus selected. And RNA-seq array data was available for all selected cases.
Raw data from EGFR-WT and EGFR-MUT cases, obtained from TCGA HT-seq files, were analyzed using Limma package of R for differential gene expression. 1598 genes were up-regulated in EGFR-MUT cases at FDR<0.05(Log2FC>1), while 663 genes were down-regulated in EGFR-MUT cases at the same threshold. EGFR was found to be up-regulated at a 2.3-fold (p = 1.69E-14) as well (Supplementary Table.2).
We further identified the biological processes enriched in the EGFR-MUT and EGFR-WT cases respectively, and enrichment analysis was employed among the differentially expressed genes using PANTHER. Genes involved in immune system process (p = 1.50E-12) and immune response (p = 4.93E-09) were enriched in EGFR-MUT LGGs, and simultaneously, genes participating in biological development were the top processes enriched in EGFR-MUT cases. Those genes related to nervous system development (p = 2.26E-22), neurogenesis (p = 4.50E-20) and generation of neurons (p = 1.30E-19) were top three processes enriched in EGFR-WT cases. Details shown in supplementary Table.3.
To meticulously analyzed the enriched gene sets, sectionalized GO analysis was performed based on the grouping strategy of dividing into Log2FC≥2 set and 2>Log2FC≥1 set. The results showed that gene sets involved in biological development were principally enriched in the Log2FC≥2 set. While immune response related gene sets were dominantly enriched in the 2>Log2FC≥1 set (Fig.2.A-B).
Gene Set Enrichment Analysis (GSEA) was also enrolled to analyze the differential expressed genes set. A pre-ranked gene list was generated based on Log2FC values from the Limma analysis, and the list were then performed GSEA pre-ranked analysis with the hallmark gene sets used as the background. Total 31 gene sets were significantly upregulated in EGFR-MUT cases and 3 gene sets were enriched in EGFR-WT cases at FDR<25%. Several immune response related gene sets such as interferon-γ response, interferon-α response, allograft rejection, TNFα signaling via NFκB, inflammatory response, IL6-JAK-STAT3 signaling, IL2-STAT5 signaling, complement and UV response were significantly enriched in EGFR-MUT cases.. While the gene sets significantly enriched in EGFR-WT cases were hedgehog pathway, pancreas β cell and β-catenin pathway. The results were in consonance with Gene Ontology enrichment analysis. Details were showed in Fig.2.C-N.
Enrichment of immune response gene sets correlate with infiltration of specific immune cell types in EGFR mutant LGGs.
As immune response related gene sets were characteristically enriched in EGFR-MUT cases, we investigated whether this enrichment trait was correlated with increased tumor infiltrated immune cells (TIIC). The immune cell infiltration levels were compared using data available at TIMER and Wilcoxon test between EGFR-MUT and EGFR-WT subset.
Consistent with the GO and GSEA results, the present results revealed a higher infiltration level of CD4+ T cells, neutrophils, macrophages, and dendritic cells in EGFR-MUT cases(Fig.3.A-D) To further explore the relationship between EGFR mutation and TIIC, the infiltration levels were analyzed according to EGFR copy number variations, by using data available at TIMER. Infiltration levels of subsets of CD4+ T cells, neutrophils, macrophages, and dendritic cells were significantly increased with EGFR level amplification cases compared to diploid cases (Fig.3.E-H). For neutrophils and macrophages subsets, infiltration levels in EGFR amplification (Copy-number > 2) cases were significantly higher than those in EGFR level-gain (2>Copy-number >0) cases. Additionally, EGFR copy-number and mRNA level were observed significantly higher in EGFR-MUT LGGs (Fig.3.I-J).
On the other hand, relationships were analyzed between EGFR mutation and genes from enriched immune response pathways. To extracted genes which actively participate in the immune response changes triggered by EGFR mutation, we chose the genes that were significantly associated with clinical outcomes both in EGFR-MUT and EGFR-WT cases(FDR<0.05). A total of 34 genes met this criterion (supplementary Table.4). Then we unexpectedly found that PD-L1 and PD-1 were elevated in EGFR-MUT LGGs at mRNA level (Fig.3.K-L). Protein-protein interaction (PPI) network analysis demonstrated that CD274 (PD-L1) was one of the hub genes among this dataset (Fig.3.M). And a Pearson correlation analysis revealed that PD-L1 expression was correlated with higher TIICs infiltrations in LGGs to a certain extent (Fig.3.N-U). A hypothesis was consequently emerged that: The increasement of both EGFR copy number and mRNA level might cause a series of biological changes including a PD-L1 elevation and rise the TIICs level, eventually change the immune microenvironment. Similar immune signature was previously reported in EGFR mutation-positive non-small-cell lung cancer (NSCLC)[28-30], in which, elevated PD-L1 in EGFR mutant tumor was consider to be an indicator of better overall response rate to PD-1 inhibitors.
Increasing infiltration of immune cells correlates with poor overall survival.
High levels of immune cell infiltration were generally considered correlated with better survival in diverse tumors. We inspected whether there was a resembling effect in the survival of LGG patients with or without EGFR mutation. Survival analysis was performed based on the grouping of EGFR mutant and infiltration levels of different types of immune cells respectively. We found a significant difference in the survival of patients with and without EGFR mutation (p <0.0001), displaying an inferior prognosis associated with EGFR mutation. And identified infiltration of CD4+ T cells, neutrophils, macrophages, and dendritic cells possess excellent biomarker potential for monitoring prognosis (Fig.4.A-D). Therefore, we particularly investigated the effects of infiltration of immune cells on clinical outcomes in all LGG cases enrolled, as well as in EGFR-MUT and EGFR-WT respectively. In the whole cohort, higher CD4+ T cells, neutrophil cells and dendritic cells infiltration levels were correlated with worse survival outcomes. Similar conclusions were also acquired in EGFR-WT cases. As for EGFR-MUT cases, macrophage cells infiltration was observed correlated with poor prognosis (Fig.4.E-H).
To further explore the possible explanation of correlation of EGFR mutations and dismal prognosis, we established COX proportional hazards models using genes enriched in EGFR-MUT LGGs or EGFR-WT LGGs and infiltration levels of all types of immune cells(details shown in Supplementary Table.5). Four genes, enriched in EGFR-MUT cases: PLSCR1, TNFAIP6, IGFBP2 and PLAT, rose the hazard ratio (HR) in EGFR-WT cases at FDR<0.05. SOCS2 on the other hand reduced the hazard ratio in EGFR-WT LGGs with a p value of 0.027. PLSCR1, a regulator of cell proliferation, maturation, apoptosis and differentiation in leukemia[31]. TNFAIP6 is constitutively expressed in adult CNS and involved in astrocyte-mediated glial formation [32]. IGFBP2 was previously reported to augment the nuclear accumulation of EGFR and to potentiate STAT3 transactivation via nuclear EGFR signaling pathway in glioblastoma[33]. PLAT, also known as tissue-type plasminogen activator (tPA), is thought to regulate vascular-genesis in tumors[34, 35]. And SOCS2, identified as a STAT3 suppressor, acts in a negative feedback loop as regulators of cytokine-triggered cell signaling[36]. Afterwards, we performed log-rank test of these covariates in COX model specifically in EGFR-MUT and EGFR-WT cases. Results showed that higher expression of IGFBP2, PLSCR1 and PLAT is correlated to poor overall survival in EGFR-WT cases, and similar effects were observed in TNFAIP6 or SOCS2 high expression cases at the window of approximately 12 years. However, opposite effects were observed for EGFR-MUT cases, highly expressed SOCS2 or PLSCR1 significantly reduced the survival risks at FDR<0.05 and associated with better outcomes (Fig.5.A-H). CELSR1, a core gene from planar cell polarity (PCP) signaling pathway regulating development of the cerebral cortex[37], which was found enriched in EGFR-WT cases and could reduce the hazard ratio (HR=0.632, p value=0.05) significantly in EGFR-MUT cases. Log-rank test showed that increased expression of CELSR1 was associated with better OS in EGFR-MUT cases, but higher level of CELSR1 indicated worse survival in EGFR-WT cases at the threshold of approximately 14 years(Fig.5.I-J).