Analysis and Construction of ceRNA Networks Reveal 4 mRNAs as Potential Biomarkers of Temozolomide-resistant Glioblastomas

Background: Glioblastoma multiform (GBM) is the most malignant tumor of central nervous system. Although temozolomide (TMZ) was reported to improve the prognosis of GBM since applied to clinical practice, its chemoresistance has become a hot research focus and even been thought as the major factor of tumor recurrence. It’s urgent to illustrate the mechanism of TMZ resistance in GBM. Methods: In this study, we retrieved and analyzed expression proles of TMZ resistant and TMZ sensitive glioma cell lines from GEO datasets, and then identied the enriched pathways of the differentially expressed genes. Moreover, genes involved in the KEGG pathways closely related to GBM from the Comparative Toxicogenomics 2020 Update Database were used for survival analysis of GBM samples from TCGA database. Finally, competing endogenous RNA networks regulating the acquired genes were constructed based on the expression proles and public database. Results (cid:0) To our knowledge, we rstly identied 4 mRNAs (LAMA1, COL5A1, ITGA2 and STK11) as potential biomarkers of TMZ resistant GBMs, then constructed their competing endogenous RNA networks. Conclusions: The results would provide theoretical basis for overcoming TMZ chemoresistance and contribute to more optimal chemotherapy of GBMs.

potential biomarkers of TMZ resistant GBMs, then constructed their competing endogenous RNA networks.
Conclusions: The results would provide theoretical basis for overcoming TMZ chemoresistance and contribute to more optimal chemotherapy of GBMs.

Background
Glioblastoma multiform is the most lethal of the primary brain tumors in adults and its prognosis is still dismal despite of the improvements of the adjuvant radiotherapy and chemotherapy in last decades [1].
Temozolomide is an oral antitumor drug approved by the US Food and Drug Administration (FDA) for the treatment of glioma. As a DNA alkylating agent, TMZ is known to induce cell cycle arrest at G2/M phase and eventually lead to apoptosis [2]. The median survival of newly diagnosed GBM were improved to 14.6 months by concomitant therapy using both radiation and TMZ, whereas 12.1 months in radiotherapy alone [3]. However, TMZ resistance is a major problem in the treatment of glioblastoma since at least 50% of treated patients do not response to it. According to different mechanisms, they were divided into two subgroups of the TMZ resistance in gliomas, the intrinsic and the adapted [4]. Acquired chemoresistance in patients and glioma cell lines received TMZ therapy was referred to the adapted. While many molecular mechanisms and pathways, such as the mutation of isocitrate dehydrogenase 1 or 2, MGMT promoter methylation, DNA mismatch repair and cancer stem cells [5], have been identi ed in the intrinsic TMZ resistance, few about the adapted resistance were thoroughly revealed.
Long noncoding RNAs (LncRNAs) are de ned as transcripts of more than 200 nucleotides that are not translated into proteins, which can function as a sponge of microRNAs (miRNAs) via competing endogenous RNAs (ceRNAs), thereby rescuing the miRNA-targeted messenger RNAs (mRNAs) [6]. The complexity and diversity of potential ceRNA interactions have scaled exponentially with the identi cation of more than 10000 lncRNAs [7]. The ceRNA crosstalk among lncRNAs, miRNAs and mRNAs were important post-transcriptional regulators of gene expression in both physiology and pathology [8]. A number of ceRNAs have been identi ed to play an important role in the initiation and progression of glioma, including the stemness, survival, apoptosis, invasion, recurrence and drug resistance [9][10][11][12][13].
However, a better understanding of the ceRNA networks regulating TMZ chemoresistance in GBM cells is still needed.
In this study, we analyzed differentially expressed mRNAs, lncRNAs and miRNAs by comparation of TMZresistant and TMZ-sensitive glioma cell lines from the public dataset. Then we screened differentially expressed mRNAs (DE-mRNAs) signi cantly related to survival prognosis and constructed competing endogenous RNA network regulating them. Ultimately, we identi ed 4 mRNAs as important biomarkers of GBM cells resistant to TMZ chemotherapy. Our study may facilitate comprehensive analysis of essential molecules, aid in the identi cation of novel mechanisms underlying GBMs' chemoresistance and reveal potential therapeutic targets for treating glioblastomas. were obtained from the public GEO database(www.ncbi.nlm.nih.gov/geo). Both datasets contained 6 samples, of which 3 were TMZ-resistant U251 cell lines and the other 3 were TMZ-sensitive U251 cell lines. Then lncRNAs and mRNAs were re-annotated according to HUGO Gene Nomenclature Community (HGNC) (http://www.genenames.org/).

2.Differential expression analysis
The raw read counts of the expression pro les were normalized and analyzed in R software (Version: 3.34.0). Differentially expressed lncRNAs, mRNAs and miRNAs were identi ed by LIMMA package (https://bioconductor.org/packages/release/bioc/html/limma.html), with a signi cance threshold of an absolute log2 fold change > 0.5 and a false discovery rate (FDR) adjusted p value less than 0.05. In addition, a volcano plot by the ggplot2 package and a heatmap by pheatmap package (Pheatmap Version 1.0.8 https://cran.r-project.org/package=pheatmap) in R were used to show the differentially expressed lncRNAs, mRNAs and miRNAs.

3.Construction of protein-protein interaction network and functional enrichment analysis
The protein-protein interaction networks were constructed by String Version 10.5 (https://string-db.org/) based on DE-mRNAs and visually displayed by Cytoscape Version 3.7.2 (http://www.cytoscape.org/).
The functional enrichment analysis of the DE-mRNAs above was performed by using the online tool DAVID (Version 6.8 https://david.ncifcrf.gov/). The terms of GO and KEGG pathway with the number of enriched counts ≥3 and p values < 0.05 were considered signi cant.

4.Screening of prognosis related DE-mRNAs
The RNA sequence data of 173 GBM samples derived from Illumina Hiseq 2000 RNA Seqencing platform were retrieved from TCGA repository(https://gdc-portal.nci.nih.gov/). Patients who had received TMZ chemotherapy were selected and the single factor cox regression analysis by R package Survival Version 2.41-1 (http://bioconductor.org/packages/survivalr/)was conducted to con rm the DE-mRNAs signi cantly related to overall survival (p value <0.05).

6.Screening of important markers of GBMs resistant to TMZ
Functional enrichment analysis of DE-mRNAs affecting prognosis signi cantly were conducted. GBM associated KEGG pathways were retrieved from the Comparative Toxicogenomics 2020 Update Database(http://ctd.mdibl.org/). DE-mRNAs related to the overlapping KEGG pathways were ltered as important biomarkers regulating GBM chemoresistance to TMZ.

2.Protein-protein interactions and the function analysis of DE mRNAs
To better identify critical genes and understand the biological function of the DE mRNAs, we established a protein-protein interaction network in STRING software (Fig 2). With the interaction scores > 0.9, 617 DE mRNAs were selected to further functional enrichment analysis by the online tool DAVID. A total of 69 biological processes and 26 KEGG pathways were enriched, of which the top 15 items, according to the p value and gene number, were shown in Fig 3 respectively. We observed the most signi cant biological processes de ned by GO were the cellular response to retinoic acid (GO: 0071300), extracellular matrix organization (GO: 0022617) and canonical Wnt signaling pathway (GO: 0060070). Simultaneously, the most signi cant KEGG pathways were the pathways in cancer (KEGG: hsa05200), focal adhesion (KEGG: hsa04510) and alcoholism (KEGG: hsa05034).
3.DE-mRNAs signi cantly affecting survival prognosis and ceRNA networks regulating them To determine the DE-mRNAs signi cantly affecting the survival prognosis of GBMs, we selected and analyzed the expression pro les and clinical information of 67 GBM patients who had received TMZ chemotherapy from the public TCGA database. Single factor cox regression analysis was carried out and the p value < 0.05 was set as selection criteria. 206 DE mRNAs were reserved for the construction of ceRNA network. As miRNA sponges, lncRNAs can isolate and bind miRNAs to regulate mRNA expression.
Differentially expressed lncRNAs and miRNAs were ltered by the predictive databases of DIANA-LncBasev2 and starBase Version 2.0, then the lncRNA-miRNA-mRNAs interactive network were constructed and displayed in gure 4. A total of 83 and 356 interaction paired relationships of lncRNA-miRNA and miRNA-mRNA were identi ed and shown respectively. In the ceRNA network, 44 DE-mRNAs were downregulated by 14 lncRNAs via 19 miRNAs, while 29 DE-mRNAs were upregulated by 29 lncRNAs via 3 miRNAs (Fig 4).

4.Screening of TMZ resistant biomarkers in GBM cells
Biological functions of the 73 DE-mRNAs (44 downregulated, 29 upregulated) acquired last step were further investigated by function enrichment analysis. A total of 10 GO biology processes and 8 KEGG pathways were obtained (Table 1). We further searched for important KEGG pathways of glioblastoma in Comparative Toxicogenomics 2020 Update Database (CTD) ( Table S1). Then we selected the 6 overlapped KEGG pathways and used the pathway-related mRNAs for survival prognosis analysis. The mRNAs expression pro les of 67 patients who had received TMZ chemotherapy (Table S2-

Discussion
Glioblastoma multiform is a grade IV brain tumor characterized by a heterogenous population of cells that are highly in ltratively, angiogenic and resistant to chemotherapy [14]. As the most widely used chemotherapy drugs, TMZ has been proved to prolong the median survival of GBM patients obviously. However, long-term survival for GBM patients remains uncommon as cells with intrinsic or adapted resistance to TMZ repopulate the tumor. Many mechanisms, such as DNA damage response pathways [15], cancer stem cells[16], microenvironment-mediated chemotherapy resistance [17], tumorderived endothelial cells[18], autophagy [19], were proved to play an important role in TMZ chemoresistance. However, the role of ceRNA networks in TMZ chemoresistance of glioblastomas was still not revealed completely.
The ceRNA hypothesis proposed in 2011 described an intricate post-transcriptional regulatory network that mainly includes lncRNAs, miRNAs, mRNAs, circular RNAs (circRNAs), and other types of RNAs [20]. A number of ceRNAs have been identi ed to play an important role in the initiation and progression of glioma. The MIR155HG/has-miR-129-5p/C1S axis promoted the proliferation, migration and invasion ability of GBM cells and was a potential marker and therapeutic target for GBM patients [21]. The circ-0043278/miRNA-638/ Homeobox A9 (HOXA9) axis had a vital function in promoting GBM progression [22]. In this study, TMZ chemoresistance related lncRNAs, miRNAs and mRNAs were screened from public databases and the crosstalk networks were predicted and constructed, in nally, the survival analysis was performed to identify the TMZ chemotherapy associated biomarkers in glioblastomas.
In general, 4 mRNAs (LAMA1, COL5A1, ITGA2 and STK11) and 3 related KEGG pathways (hsa04510: Focal adhesion, hsa04512: ECM-receptor interaction, hsa04151: PI3K-Akt signaling pathway) were con rmed to play an important role in TMZ chemoresistance. LAMA1 and COL5A1 were upregulated, while ITGA2 and STK11 were downregulated in TMZ-resistant glioma cells. The gene expression of laminin subunit alpha 1 (LAMA1) was signi cantly higher in glioblastomas than in non-neoplastic white matter, veri ed by immunohistochemistry and real time quantitative PCR [23]. By analyzing the expression pro les of recurrent low grade glioma (LGG) and primary LGG, collagen type V alpha 1 chain (COL5A1) was identi ed as one of the hub DEGs and might participate in the pathological process of recurrence [24]. Furthermore, COL5A1 was a metastasis-associated gene in clear cell renal cell carcinoma and immune cell in ltration related gene in head and neck squamous cell carcinoma or gastric cancer [25][26][27]. Integrin subunit alpha 2 (ITGA2) was identi ed as a novel molecular target of GBM and its antibody blockage signi cantly impeded GBM cell migration [28]. Serine/threonine kinase 11 (STK11) is one member of the serine/threonine kinase family, which is involved in regulating cell polarity, apoptosis, and DNA damage repair [29].
In addition to mRNAs, we screened 34 lncRNAs and 21 miRNAs through the ceRNA network. As shown in gure 5, miR-9-5p is one of the hub miRNAs in the network and may sensitize GBM cells to TMZ by downregulation of LAMA1 and COL5A1. Consistent with our prediction, miR-9-5p was reported to inhibit GBM cells proliferation by downregulating forkhead box P2 (FOXP2) [30]. Of the 34 lncRNAs, LUCAT1, TCL6, MIR4458HG and LINC00114 may play a vital role in TMZ chemotherapy resistance of GBM cells. Yansheng gao et al. revealed that LUCAT1 could act as a oncogenic lncRNA and promote glioma tumorigenesis via regulating miR-375 [31].
Although we screened 4 genes as important biomarkers of TMZ chemoresistance in GBM and established the regulatory ceRNA network. There are a few limitations to this study. First, veri cation of ceRNA axis based only on public database is insu cient, further cell experiments are needed to con rm our results. Second, we evaluated genes' effects on TMZ chemoresistance only by survival prognosis analysis based on TCGA database, prospective researches are needed to validate our results. Third, 6 KEGG pathways were overlapped between function enrichment analysis of prognosis associated DE-mRNAs and CTD database. The TMZ resistant cell lines' expression levels of genes involved in 3 pathways (hsa04540: Gap junction, hsa04912: GnRH signaling pathway, hsa05202: Transcriptional misregulation in cancer) contradicted their prognosis in GBM patients. So we selected 3 KEGG pathways for further analysis.

Conclusions
To our knowledge, we rstly identi ed 4 mRNAs as important biomarkers of TMZ resistant GBM, then constructed their competing endogenous RNA networks. Although further experiments are needed to con rm our results, this study would provide theoretical basis for overcoming TMZ chemoresistance and contribute to more optimal therapy of GBM.

Availability of data and materials
The datasets used and/or analysed during the current study are available at https://www.ncbi.nlm.nih.gov/gds/ (GSE100736 and GSE100775) or from the corresponding author on reasonable request.     The differentially expressed lncRNAs, miRNAs and mRNAs regulatory networks. Square, triangle and circle represent lncRNA, miRNA and mRNA, respectively, and the color depicting the expression level in TMZ-resistant U251 cell lines. CeRNA regulatory network and KEGG pathways of the 4 biomarkers in TMZ resistant GBMs. Square, triangle, circle and rhombus represent lncRNA, miRNA, mRNA and KEGG pathway, respectively. The color depicts the expression level of differentially expressed RNAs.

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