Identication Of Immune Networks Inuenced By Mir-200c In IDH1-Mutant Lower-Grade Gliomas Using WGCNA

Gliomas are the most common intracranial tumors with a highly invasive nature and poor prognosis. Different stages of gliomas have distinct alterations in genome, genetic and epigenetic changes. Here, we performed enrichment analysis and weighted gene correlation network analysis (WGCNA) on the Cancer Genome Atlas (TCGA) RNA-seq dataset to identify core genes and potential pathways involved in Isocitrate Dehydrogenase 1 (IDH1) mutant lower-grade gliomas (LGG). The IDH1-mutant LGG patients were found to frequently co-occur with abundant genetic alterations, surprisingly, the miR-200c/141 remain unchanged. Although located in close proximity, the expression level of miR-200c was more affected by methylation level and copy number variation than that of miR-141 in IDH1-mutant LGG. GO analysis indicated the high correlation between miR-200c and immune response. Gene set enrichment analysis (GSEA) suggested that leukocyte activation, inammatory response, lymphocyte activation and adaptive immune response were signicantly up-regulated in miR-200c high-expression group. Enrichment analysis demonstrated that leukocyte activation, immune response, lymphocyte proliferation, cell-cell adhesion and other items were participants in IDH1-mutant LGG samples. Furthermore, WGCNA yielded 33 signicant modules, including the yellow module which contains genes involved in the immune response. Taken together, these results proposed a promising new therapeutic strategy for treating LGG gliomas through a miR-200c-mediated immune response pathway.


Introduction
Glioma, which originates from the glial cells, is the most common and aggressive primary brain tumor, accounting for the majority of brain-related deaths [1,2]. Based on the 2007 World Health Organization (WHO) classi cation of central nervous system (CNS) tumors, gliomas were categorized by cell type and graded from I to IV (lowest to highest invasiveness) [3]. Benign Grade I glioma, as well as grade II glioma with a survival time of more than 5 years, were classi ed as low-grade gliomas. While grade III glioma with a survival time of 2-3 years and IV gliomas (Glioblastoma, GBM) were de ned as high-grade gliomas [4]. Traditionally, grade II gliomas were considered to have a better prognosis than grade III gliomas [5]. However, with the identi cation of key molecular alterations such as isocitrate dehydrogenase (IDH) 1/2 mutations in more than half of the grade II and III gliomas, WHO introduced a new classi cation of the CNS tumors in 2016, and the term "lower-grade gliomas (LGG)" that covers grade II and III gliomas gradually replaced the term "low-grade gliomas" [6]. IDH1 mutations were rst discovered in 12% of GBM patients through an integrated genomic analysis of human glioblastoma multiforme in 2008 [7].
Subsequently, IDH1-mutations were identi ed in over 70% of grade II/III glioma and secondary GBM cases [8]. Despite considerable efforts in the past decades toward treating gliomas, the survival rate of patients remains dismally low. There are hundreds of molecular alterations in glioma. Although previous studies have shown that IDH1 mutation is an independent favorable prognostic marker in gliomas [9], our knowledge about its initiation and progression is incomplete. The pharmacologic manipulation of key molecular pathways to improve patients' survival rate remains to be further explored.
MicroRNAs (miRNAs) are small non-coding RNA molecules that have a signi cant regulatory capacity on the expression of downstream genes [10,11]. MiR-200 is a family of tumor suppressor miRNAs consisting of ve members (miR-200a, miR-200b, miR-200c, miR-429, and miR-141), which are signi cantly involved in the inhibition of epithelial-to-mesenchymal transition (EMT), repression of cancer stem cells selfrenewal and differentiation, modulation of cell division and apoptosis, and reversal of chemoresistance. The targets of the miR-200 family include the 124kDa transcription factors zinc nger E box-binding homeobox 1/2 (ZEB1/2), which act as markers of EMT and increase metastasis by repressing E-cadherin and Lgl-2 [12]. Among the miR-200 family members, miR-200c is of particular interest for human health and disease, and it has recently emerged as an important regulator of tumorigenicity and cancer metastasis. However, its role in IDH1-mutant LGG is not de nitively clear.
In this article, we performed enrichment analysis and weighted gene correlation network analysis (WGCNA) to explore the mechanisms of miR-200c in IDH1-mutant samples and compared its prognostic value in LGG. The results showed a high correlation between miR-200c and immune response including leukocyte activation, in ammatory response, lymphocyte activation, and adaptive immune response in IDH1-mutant LGG samples, which give us a hint to treat LGG using immune checkpoint receptors. Additionally, WGCNA yielded yellow module is highly associated with miR-200c/miR-141 and contains genes involved in the above immune response. Taken together, these results proposed a promising therapeutic strategy for treating LGG through targeting a miR-200c-mediated immune response pathway.

Data acquisition and processing
The RNA-seq datasets, DNA methylation datasets, and clinical datasets from patients with LGG were acquired from the UCSC Xena browser (https://xenabrowser.net/) [13]. The expression of miR-200c/141 in IDH1 mutant LGG and the corresponding normal tissues were obtained from The Cancer Genome Atlas (TCGA) datasets and Genotype-Tissue Expression Portal [14]. Ensembl database (http://asia.ensembl.org/index.html, GRCH38.p13) was applied to integrate ID into o cial gene symbols. R Maftools package was used to analysis of somatic variants in IDH1 mutant. Kaplan-Meier method was used for survival analyses, with a P-value generated via log-rank test.

Differential Expressed Genes (DEGs) Visualization
The DEGs analysis of samples with the highest expression of miR-200c (n=50) and the lowest expression (n=50) was identi ed using the R limma package, with a Q value < 0.05 and Abs (FC) > 2 as cutoff criteria. Volcano plots were generated in R software using the VolcanoPlot package. The Heatmap.2 package was used to visualize the expression heatmap of DEGs. Gene Ontology analysis of DEGs was performed using Panther GO (http://geneontology.org/page/go-enrichment-analysis). The ggplot2 package was applied to visualize the results. GSEA analysis of the highest and the lowest expression of miR-200c was performed by using GSEA 3.0 (http://www.broadinstitute.org/gsea/). The heatmap of different GO clusters was generated by Heatmap.2 package.
2.3 WGCNA of RNA-Seq Data from LGG Patients WGCNA R package was performed to, rather than only nd the DEGs, identify the associations between miR-200c/141 expression and each module by the 'p. weighted' function. A weighted network adjacency was de ned by increasing the parameter β which is used to optimize the correlations between genes. The power of β = 4 and scale-free R 2 = 0.95 were applied as the soft-thresholding parameter. A hierarchical clustering tree was constructed, with different branches of the tree representing different gene modules. The adjacency matrix was transformed into a topological overlap matrix (TOM). Genes were divided into different gene modules based on the TOM-based dissimilarity measure. Then, the module-trait relationship between the clinical data (gender, laterality, neoplasm histologic grade, histological type), the miR-141/200c, and different modules were analyzed. The gene interaction network was constructed by the String database (https://string-db.org/), and the interactive relationship of all the yellow model networks was calculated and visualized by Cytoscape 3.4.0. GO analysis of the genes in the yellow module was performed using the clusterPro ler package and Metascape [15,16].
The locations of miR-200c/141 determine the target speci cation. Among the targets, miR-200c/141 were reported to target the ZEB1/2 transcript factor for degradation [17,18]. To demonstrate whether the expression of ZEB1/2 impacts miR-200c/141 in IDH1-mutant tumors, we used TCGA datasets to analyze the expression pattern of miR-200c/141and ZEB1/2. Results showed that compared with IDH1 wild-type tumors, the expression of miR-200c/141 were both down-regulated, while the expression of ZEB1/2 was up-regulated in IDH1-mutant tumors, suggesting that ZEB1/2 might be the functional transcriptional suppressor of miR-200c/141 (Figure 1c-f). This inverse relationship between ZEB1/2 and miR-200c/141 has been previously identi ed from a panel of human cancer cell lines [18]. Taken together, these results suggested that miR-200c and miR-141 play important roles in IDH1-mutant glioma patients' OS via a double-negative feedback loop with ZEB1/2.

MiR-200c is preferentially up-regulated in IDH1-mutant
LGG samples  and miR-141 are both located on chromosome 12p13 within the same microRNA cluster according to chromosome locations [19]. We performed Pearson Correlation Coe cient analysis to identify the association between miR-200c and miR-141, and the results showed a strong positive correlation between the expression level of miR-200c and miR-141 (Figure 2a, b). Previous studies of DNA promoter alterations from TCGA revealed that the miR-200c/141 cluster is epigenetically regulated by DNA methylation [20]. To determine the mechanism of miR-200c/141 in regulating IDH1-mutant LGG, we investigated the relationship between DNA methylation and miRNAs expression. For this purpose, we extracted a 2kb sequence from the upstream region of miR-200c and miR-141 to analyze the DNA methylation status of both genes. Six CpG sites were identi ed using UCSC Genome Browser on the promoter region (Figure 2c). Since one probe data was inaccessible, we plotted a heatmap relating the other 5 CpG sites (cg16642299, cg24702147, cg00366413, cg27534624, and cg22413603). All probes  (Figure 3c). These differentially expressed genes were analyzed by GO analysis. The biological process showed that miR-200c related genes are of the highest association with neutrophil activation involved immune response. The other processes include neutrophil activation, T cell activation, neutrophil degranulation, neutrophil-mediated immunity, regulation of lymphocyte activation, suggesting a strong correlation between miR-200c and immune response (Figure 4a, b) GSEA indicated that leukocyte activation, in ammatory response, lymphocyte activation, and adaptive immune response were signi cantly up-regulated in miR-200c high-expressed group, suggesting the role of miR-200c in activating immune response to prevent diseases (Figure 4c).

WGCNA on RNAseq in IDH1-mutant LGG samples
WGCNA is a data mining method that has been widely used to study biological co-expressed gene networks by non-scale weight network construction and module detection [21]. To verify the gene network that was highly associated with miR-200c expression, WGCNA was applied to further analyze IDH1mutant LGG RNAseq data. We rst downloaded RNAseq data from the TCGA database, and the samples of IDH1 mutant were clustered using the average linkage method and Pearson's correlation method. The co-expression analysis was carried out to construct the co-expression network. In this study, the power of β = 5 (scale-free R 2 = 0.8) was selected as the soft-thresholding parameter to ensure a scale-free network.
The dendrogram shown in Figure 5a (Figure 5d). The results suggested that genes that highly correlated with the yellow module were also strongly associated with miR-200c expression.

Enrichment analysis of genes highly associated with miR-200c expression in IDH1-mutant LGG samples
The network connections in the yellow module, with weight parameters higher than 0.19 from WGCNA, were loaded into Metascape and analyzed. A variety of proteins were included in the network, such as CD180, ITGAL, SYK, ITGB2, CD53, and LAPTM5 (Figure 6a). High inner connectivity of the gene cluster was depicted by the edge sizes. Finally, the genes from the yellow module with a coe cient higher than 0.5 were selected and underwent GO enrichment analysis (Figure 6b). Leukocyte activation, immune response, lymphocyte proliferation, cell-cell adhesion, and other items were enriched in the analysis, indicating the key role of miR-200c in immune responses in IDH1-mutant LGG tissues (Figure 6b). Through Pearson correlation, miR-200c was found signi cantly correlated with molecules associated with its function including miR141, MLPH, TTC39A, C1orf64, and GSTM3, and other immune checkpoint molecules including miR141, CTLA4, LAG3, CD160, and CD276 (Figure 6b).

Discussion
Glioma is a common nervous system cancer that represents di culties in completely surgical removal [22]. Along with the resistance to radiotherapy and chemotherapy, glioma has become one of the most malignant tumors [23,24].
LGGs are in ltrative gliomas that most often occur in the adult cerebral hemispheres, a subset of LGGs stays stable for several years, whereas some progress into GBM within a few months, causing a dismal prognosis. Hence, understanding the molecular mechanisms of glioma initiation and progression is important in developing effective therapeutic strategies for LGG patients. In this study, we identi ed miR-200c, which belongs to the tumor suppressor miR-200 family, is crucial for IDH1-mutant LGG patients through activating an immune response.
MiR-200c belongs to the miR-200 family [25]. In recent years, an increasing number of studies have shown that miR-200c expression inhibits the progression, invasion, EMT and chemoresistance of a variety of tumors, including non-small cell lung cancer, breast cancer stem cells, prostate cancer and ovarian cancer [26,27] [28, 29]. Furthermore, miR-200c/141 were reported to inhibit the expression of transcription factors ZEB1 and ZEB2, which act as markers of EMT, resulting in the downregulation of E-cadherin expression and reduction of cell adhesion and nally promoting the invasion and metastasis of tumor cells [30,31]. These results indicated that despite the mutually genetic alterations with IDH1 mutation existing in LGG patients, no miR-200c mutation was detected, which eventually resulted in better overall survival. It is consistent with our results, that is, compared with IDH1 wild-type tumors, IDH1-mutant tumors with a relatively higher ZEB1/2 expression showed a lower expression of miR-200c/141, indicating the vital role of miR-200c as a tumor suppressor through regulating the expression of ZEB1/2. The opposite expression pattern of miR-200c/141 and ZEB1/2 demonstrates the negative regulatory feedback loop between miR-200 family members and ZEB TFs, future studies will be needed to test in LGG samples. IDH1 mutations are closely related to epigenetic variations. Epigenetic mechanisms, especially DNA methylation, play a decisive role during tumorigenesis and EMT and contribute to the regulation of key factors involved in this process. The hypermethylation in regulatory regions of many tumor suppressor genes leads to their transcriptional silencing (e.g. E-cadherin [32]). MiR-200c and miR-141 are both located on chromosome 12p13, however, Niu et al. identi ed a transcriptional unit containing miR-200c, but not miR-141, within the CpG Island at position -285. The result of an alternative splicing variant exists to splice out miR-141 and conserve the miR-200c transcript was supported by an in-silico inspection, upon which a high-score splice donor site (boundary exon/intron) was detected between the two miRNAs in position +373 [20]. This data might explain the observation previously made by us and others [20,21], that miR-200c is usually expressed at higher levels and is more affected by methylation level and copy number variation than miR-141, despite their close location.
Hypermethylation of CpG islands in promoter regions contributes to the miR-200c expression [20]. According to the epigenetic methylation level, IDH1-mutant LGGs are normally classi ed as methylationrich. In our case, the methylation level of the miR-200c gene was detected by the six probes on the CpG island before the transcription start site of miR-200c. Moreover, there is a signi cant association between methylation level and ampli cation-type alteration, which were both associated with the expression level of miR-200c. To further understand the molecular mechanism behind the effect of miR-200c in IDH1mutant LGG tissues, we performed RNA-seq analysis, downstream GO analysis, and GSEA. The results illustrated that miR-200c might be an important regulator of the immune response, which was consistent with the previous study that miR-200c can control the in ammatory process and the response to an infection [33][34][35].
It was shown that miR-200c has the ability to modify the e ciency of TLR4 signaling through the MyD88dependent pathway, thus in uencing host innate defenses against microbial pathogens [35]. Rogers, T.J., et al. illustrated that miR-200c led to reduced production of the immune-suppressive metabolite kynurenine directly by targeting TDO2. Furthermore, in addition to reversing a classical EMT signature, miR-200c represses many genes encoding immune-suppressive factors containing PD-L1/2, HMOX1, and GDF15 [34]. It is also recognized that in ammatory-related miR-200c can sustain innate immune responses through LPS-dependent mechanisms [33], strengthening our result that miR-200c plays a vital role in immune responses, providing a therapeutic strategy via a miR-200c-modulated immune response.
These results were further con rmed by WGCNA and downstream enrichment analysis, showing that the genes in the yellow module played essential roles in immune responses. Besides, the genes highly correlated with miR-200c expression showed intensive connection by the visualization of metascape. A variety of the hub-genes were immune proteins, e.g., CD180, ITGAL, SYK, ITGB2, CD53, and LAPTM5, indicating the high association between miR-200c and immune systems.

Conclusions
In conclusion, all the data presented above suggested that miR-200c played an essential role in the physiological process of IDH1-mutant LGG. Based on the molecular and clinical analysis, our study elucidated that miR-200c linked tightly to immune responses in IDH1-mutant LGG tissues, and was involved in the survival of IDH1-mutant LGG patients through targeting intermediates in signaling pathways essential for innate immune responses.

Availability of data and materials
The datasets generated and analyzed during the current study are available in the TCGA Research Network: https://www.cancer.gov/tcga. All data generated or analyzed during this study are included in this published article.