Identication of ncRNA-mediated prognostic Value and Immune Inltration of MTHFD2 in Bladder Cancer

Background: Methylenetetrahydrofolate dehydrogenase 2 (MTHFD2), one of mitochondrial enzymes, is involved in folate and nucleic acid metabolism and maintains the cellular redox balance. However, the function of MTHFD2 in bladder cancer is still poorly characterized. This study was designed to elucidate the effect and regulatory mechanism of MTHFD2 on bladder cancer cells and explore the relationships between MTHFD2 and immune cell inltration in tumor microenvironment (TME). Methods: The data from Oncomine, TIMER, The Cancer Genome Atlas (TCGA) and The Human Protein Atlas (HPA) database were extracted to evaluate the expression of MTHFD2 and its prognostic role in pan-cancer, especially in bladder cancer. Enrichment analyses, were utilized to explore the underlying cellular mechanisms. The ncRNA regulatory axis was established by Starbase database, and the PPI network was constructed by Cytoscape software. Ultimately, the relations between the expression of MTHFD2 and immune cell inltration in bladder cancer was indicated by TCGA and TIMER databases. Results: Our results demonstrated that MTHFD2 expression was generally up-regulated in pan-cancers and its high expression was correlated with poor prognosis of patients with bladder cancer. Specically, our study revealed that MTHFD2 was a powerful risk factor and involved in the tumor development of bladder cancer. Furthermore, hsa_circ_0046140 and hsa_circ_0006769/miR-383-5p/MTHFD2 axis could also play a signicant role in tumorigenesis. Ultimately, a strong correlation was observed between MTHFD2 expression and various immune cell inltration, which implied that MTHFD2 might serve as an agent in tumor immunotherapy. Conclusion: MTHFD2 is widely overexpressed in pan-cancers, and its expression is related with the prognosis of patients and the multiple immune cell inltrates in TME. Besides, hsa_circ_0046140 and hsa_circ_0006769/miR-383-5p/MTHFD2 axis are implicated

MTHFD2 as a mitochondrial enzyme implicated in folic acid metabolism [7], could catalyze dehydrogenation of CH2-THF and the generation of CHO-THF, which was mediated by NADH/NAD+ in the process of folate metabolism [8]. The one carbon unit produced in this metabolic process can be utilized to synthesize purine and thymidylic acid [9,10], providing the key raw materials for intracellular nucleic acid synthesis. Meanwhile, as a NADP+ dependent enzyme, MTHFD2 can also change the level of cytoplasmic NADPH/NADH, thus affecting the redox balance of cells [11]. Down-regulation MTHFD2 could reduce the level of intracellular NADPH/NADP+ and glutathione, which suppresses the proliferation and invasion of tumor cells [12,13]. Cumulative evidences have indicated an association between MTHFD2 and several tumor progression, involving lung cancer [14], breast cancer [15], head and neck cancer [16] and kidney cancer [17]. MTHFD2 is principally expressed in mitochondria of embryonic cells, and could not be detected in most healthy adult tissues [18,19], which provides a rationale for targeting MTHFD2 in cancer therapy.
Non-coding RNAs (ncRNA), including miRNAs, lncRNAs, and circular RNAs (circRNAs), are considered as special types of RNA. Although ncRNA is incapable of encoding or translating into proteins, it ubiquitously exists in the living organisms, regulate gene expression and contribute to disease pathogenesis [20]. For example, ncRNAs could trigger gene silencing by mRNA degradation and remodel the chromatin structure, leading to the failure of protein translation. Furthermore, ncRNA can also be involved in cis and trans gene regulation, and thereby mediating the up or down-regulation gene expression [21]. Therefore, ncRNA could exert regulatory role in the proliferation, differentiation, apoptosis and so on [22]. However, the potential mechanisms and the ncRNA-mediated expression regulation network of MTHFD2 in bladder cancer have not been elucidated.
The TME is implicated in tumorigenesis, progression and invasion, and therefore, it could generate a large impact on the treatment e cacy of cancer therapy [23]. The in ltrating immune cells constitute a signi cant proportion in TME [24]. The immune cell in ltration in TME is widely assumed to be associated with tumor proliferation and migration. Immunotherapy, such as anti-PD-1/PD-L1 and anti-CTLA4 inhibitors, functioning by activating immune cells to eliminate tumor cells, has been extensively applied for the treatment of several tumors and achieved satisfactory e cacy. Hence, it is extremely urgent to illustrate the speci c landscape of immune cells in ltration in bladder cancer, in order to improve the e cacy of immunotherapy.
The function of MTHFD2, and the relations between the its expression and immune cell in ltrates in bladder cancer were still largely unknown. Here, we used bioinformatic methods to explore the associations between its expression and the prognosis of the patients. Hereafter, we investigated the biological function of MTHFD2, and ncRNAs that might be involved in regulating MTHFD2 expression. Ultimately, we revealed the relationships between MTHFD2 and immune cell in ltration. Our ndings provided a new insight into the function of MTHFD2 and a potential curative treatment for bladder cancer.

Materials And Methods
Data source from TCGA and GEO database RNA-seq and clinical data of bladder cancer patients were downloaded from the UCSC (http://xena.ucsc.edu/public). The GEO is an authoritative gene expression public database (https://www.ncbi.nlm.nih.gov/geo/). GSE13507 dataset in the GEO database was selected to evaluate the differential expression pattern of MTHFD2 in bladder cancer.

Pan-cancer analysis
Oncomine database (https://www.oncomine.org/resource/login.html) was used to explore differential expression pattern of MTHFD2 in multiple tumor tissues [25]. The TIMER (https://cistrome.shinyapps.io/timer/) was utilized for analyzing differential expression of genes and the in ltration of immune cells in a wide range of cancers tissues [26]. The differential expression analysis of MTHFD2 in pan-cancer was conducted by using TIMER and TCGA database.
The Human Protein Atlas database HPA (https://www.proteinatlas.org/) provides information of human samples, including cell, tissue and pathology atlas [27]. The immunohistochemistry staining data of MTHFD2 was retrieved from it to veri ed the differential expression pattern of MTHFD2 between the normal and tumor groups.

Enrichment analysis
Gene Set Enrichment Analysis (GSEA) was conducted by using GSEA software

MiRNA prediction analysis
StarBase is a database which was built for exploring ncRNA-related studies [29]. It contains seven targeted gene predicting tools (PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan). We chose the miRNA for further analysis only when it was included in at least two programs. Tarbase database was one of the DIANA tools, which was used to predict the relations between miRNA and mRNA [30]. MiRPathDB database was applied to explore the pathway involved the target miRNA [31]. The differential expression analysis and correlation analysis of miRNA was conducted by using R (v.4.0.3).

CircRNA prediction analysis
StarBase database was utilized to predict the relationships between the circRNA and miRNA. Cancerspeci c circRNAs database (CSCD) was employed to constructed the predicted circRNA [32].
PPI interactive network PPI interactive network was built by Cytoscape (Version 3.8.2) software. TISIDB TISIDB database is employed for accessing the correlation between the tumor cells and immune in ltration [33]. We analyzed the associations between MTHFD2 expression level and in ltration of immune cells in multiple tumors.

Statistical analysis
All statistical analysis was conducted by R (v.4.0.3). Kaplan-Meier curve, ROC curve and Nomogram chart were established by using the survival package and survminer package, the pROC package, and rms package and survival package, respectively. Univariate and multivariate cox regression and logistic regression test were employed to determine the prognostic factor of bladder cancer by survival package.

MTHFD2 expression and its association with patients' survival in pan-cancer
We initially identi ed the expression of MTHFD2 in multiple tumors were elevated except in leukemia ( Figure 1A). To evaluate the differential expression pattern of MTHFD2 between the normal and tumor groups, we used the TIMER and TCGA pan-cancer databases to further explored the expression of MTHFD2 in pan-cancer. The ndings demonstrated that MTHFD2 was overexpressed in 18 kinds of tumors compared with corresponding normal tissues ( Figure 1B Upregulation of MTHFD2 and its clinical signi cance in bladder cancer Speci cally, we analyzed the MTHFD2 expression in bladder cancer by TCGA BLAC database. The results revealed that mRNA MTHFD2 expression was distinctly increased in bladder cancer tissues, which agreed with the prior ndings (p<0.001) ( Figure 4A-B). Furthermore, the upregulated MTHFD2 expression in bladder cancer was correlated with higher T stage, M stage, pathologic stage and histologic grade ( Figure  4C-4G). Data from GEO database also manifested that MTHFD2 expression was enhanced in bladder cancer samples in comparison to normal tissues (p<0.05) ( Figure 4H). Immunohistochemical staining data from HPA also revealed MTHFD2 protein was upregulated in bladder cancer tissue ( Figure 4I-4J). We evaluated the diagnostic value of MTHFD2 and the area under the ROC curve (AUC) was 0.741 ( Figure   4K), which further veri ed the prognostic value of MTHFD2 in bladder cancer. Hereafter, we established a nomogram model by combining MTHFD2 expression and different clinicopathological characteristics (T stage, N stage, M stage and histologic grade) to predict the OS ( Figure 4L). The higher score, the worse prognosis in the nomogram model. The ndings supported that nomogram may be a reliable model for the prediction of bladder cancer patients' survival by the performance of Calibration curve ( Figure 4M).
Univariate Cox analysis showed that TNM stage, pathologic stage, age, subtype, lymph vascular invasion and elevated MTHFD2 were the contributing factors in patients' OS (p<0.05) (Figure5A). Multivariate Cox analysis revealed that MTHFD2 was a pivotal risk factor for OS (p=0.038) ( Figure 5B). Univariate logistic regression analysis demonstrated that the poor prognosis factors included T stage, pathologic stage, radiation therapy, race and subtype were associated with higher expression of MTHFD2 (p<0.05) ( Table   1).

Functional enrichment analyses
GSEA was carried out based on mRNA expression data of TCGA BLCA samples to explain the function of MTHFD2 in bladder cancer biological process. The results revealed that elevated expression of MTHFD2 was positively associated with intracellular functions including cell cycle, DNA replication, mismatch repair and DNA repair, and might be involved in the MYC, MAPK, Wnt, NF-κB and ATM signaling pathways ( Figure 6). Taken together, the GSEA results demonstrated that MTHFD2 could contribute to bladder cancer tumorigenesis.
We subsequently selected the top 300 genes which were the most associated with MTHFD2 for GO and KEGG pathway analysis. The results revealed that MTHFD2 was closely correlated with cell proliferation, including the cell cycle regulation, cell metabolic process, DNA repair, kinase activity and DNA replication ( Figure 7A-D). We subsequently constructed a network by selecting the tumor-related gene terms from the enrichment analyses. The network exhibited multiple genes which were prominently enriched in the gene terms ( Figure 7E). Prediction and analysis of miRNA regulating MTHFD2 Cumulative evidences indicated that ncRNA could regulate gene expression level, especially the miRNA.
Firstly, we used Starbase database and Tarbase database to predict 146 and 94 miRNAs that might potentially bind to MTHFD2. Subsequently, we selected the intersection of miRNAs predicted by the two databases by using Venn diagram and screened 30 identical miRNAs ( Figure 8A). Ultimately, a regulatory network of miRNA-MTHFD2 was construct by Cytoscape software ( Figure 8B). Yellow represents miRNA up-regulated in bladder cancer, pink represents down-regulated and orange represents no statistical difference. Additionally, we analyzed the intracellular functions of the 30 screened miRNA ( Figure 8C). According to the mechanisms of miRNA regulating target gene expression, we selected the downregulated miRNA, such as miR-383-5p and miR-1-3p in bladder cancer for further analysis ( Figure 8D-8E). Correlation analysis showed that miR-383-5p was inversely associated with the MTHFD2 expression ( Figure 8F-8G). Therefore, miR-383-5p may be a crucial miRNA regulating the expression of MTHFD2.
Prediction and analysis of circRNA regulating miR-383-5p To further explore the ncRNA regulatory network, the Starbase database was employed to predict upstream circRNAs of miR-383-5p. 12 circRNA which were most likely to be implicated in the regulation of miR-383-5p were selected. In the meanwhile, a regulatory network diagram was established by Cytoscape. The structural patterns of the 12 kinds of circRNAs were exhibited in Figure 9B-9M. In order to further elucidate the cellular functions of 12 circRNAs and their roles in bladder cancer, we analyzed the differential expression pattern and survival of circRNA encoding genes in bladder cancer ( Table 2). The results demonstrated that 4 circRNAs were up-regulated in bladder cancer, and among them, hsa_circ_0046140 and hsa_circ_0006769 were associated with poorer prognosis. Hence, we assumed that hsa_circ_0046140 and hsa_circ_0006769 might be implicated in the modulation of miR-383-5p.

MTHFD2 expression and tumor in ltrating immune cells
Some immunological studies were conducted to assess whether expression of MTHFD2 could affect the immune cell in ltration. We found that higher MTHFD2 expression group exhibited a higher level of immune cells in ltration, such as B cells, CD8 T cells, cytotoxic cells, eosinophils, macrophages, neutrophils, T cells, T helper cells (Th1 cells, Th2 cells), Tem, TFH, Tgd, and Treg cells (p<0.05) ( Figure  10A). The correlation analysis indicated that MTHFD2 expression is predominantly positive correlated with immune cell in ltration in bladder cancer, including Th2 cells, Th1 cells and macrophages ( Figure  10B). TISIDB database was used to further validate the correlations between MTHFD2 expression and 28 kinds of tumor in ltrating immune cells in pan-cancer ( Figure 10C). Besides, TIMER and TCGA database were employed to further examine the underlying relations between the MTHFD2 expression and immune markers ( Table 3). The ndings manifested a strong association between the MTHFD2 expression level and immune markers of CD8+T cells, Th1 cells, M2 macrophages, dendritic cells and NK cells in bladder cancer.
MTHFD2 expression is correlated with immune checkpoints PD-1/PD-L1 and CTLA-4 were crucial immune checkpoints, which partially account for the poor e cacy of tumor therapy. Therefore, we employed TIMER and TCGA database to evaluate the association between MTHFD2 and PD-1/PD-L1 or CTLA-4. The ndings showed that MTHFD2 expression is highly positive correlated with PD-1 ( Figure 11A-11B), PD-L1 ( Figure 11C-11D) and CTLA-4 ( Figure 11E-11F) in bladder cancer. These results suggest that immune checkpoints may be implicated in the tumorigenesis and development of MTHFD2 mediated bladder cancer.

Discussion
MTHFD2 is a folic acid metabolism enzyme which could mediate mitochondrial one-carbon production [7]. The generation of one carbon unit in folic acid metabolism is a key step for the nucleotide and DNA synthesis [34]. The alteration in the metabolic process of one carbon unit will exert an effect on cellular nucleic acid replication. Cumulative evidences has proved that MTHFD2 functioned in tumor cell growth, the development of drug resistance and tumor stem cell-like characteristics by affecting folate metabolism, one carbon unit synthesis and intracellular redox balance [35]. It is reported that MTHFD2 could facilitate the methylation level of hypoxia-inducible factor (HIF)-2a mRNA and elevate its translational level, which eventually leaded to the activation of glycolysis and increase of MTHFD2 expression [17]. The positive feedback cycle could enhance metabolic reprogramming and tumor development.
The reactive oxygen species (ROS) have a cytotoxic effect on tumor cells, and therefore, the increase level of ROS in tumor cells was considered as a protective factor for cancer patients [36]. As an NADP+ dependent enzyme, MTHFD2 can regulate intracellular level of NADPH/NADH and affect intracellular redox balance [37]. By downregulating the expression of MTHFD2, the intracellular level of NAPDH and glutathione could be reduced while the intracellular oxidation level could be elevated, contributing to the elimination of tumor cells. The latest study found that MTHFD2 could modulate the stabilization of ubiquinol-cytochrome C reductase core protein 2 (UQCRC2), which belongs to the Complex III family. Moreover, MTHFD2 exhaustion could result in mitochondrial functional disturbance by restraining the activity of Complex III via regulating UQCRC2 expression [38]. However, the function of MTHFD2 in modulating the proliferation of bladder cancer cells through altering the level of ROS and its potential mechanisms is still required to be investigated.
In the present study, we performed a bioinformatics analysis of RNA sequencing data retrieved from the TCGA and GEO databases. Firstly, we found that MTHFD2 was generally overexpressed (Figure 1), and higher MTHFD2 expression was related with worse overall survival in in pan-cancers (Figure 2-3). Especially, we focused on the expression level of MTHFD2 in bladder cancer and its effect on patient's prognosis. The ndings demonstrated that MTHFD2 expression was elevated in bladder cancer ( Figure  3A-3B) and overexpression of MTHFD2 is associated with unfavorable clinical pathological features, shorter survival time as well as poorer prognosis ( Figure 3C-3H). In addition, we further explored the function of MTHFD2 by enrichment analysis in bladder cancer. GSEA enrichment test manifested that the overexpression of MTHFD2 was related with biological processes including cell cycle, DNA replication and repair, and the activation of signaling pathway including MYC, MAPK, Wnt, NF-κB and ATM ( Figure  6). GO and KEGG analysis revealed that MTHFD2 might be involved in the process of tumor initiation, which was consistent with the previous results (Figure 7). The study implied that MTHFD2 might be a possible prognostic biomarker and promising targets for the treatment of bladder cancer.
It is con rmed that ncRNA modulated the gene expression through the ceRNA mechanism, which plays an important role in tumor initiation [39]. To explore miRNA that possibly regulate MTHFD2 expression, we used Starbase databases and Tarbase databases to predict miRNA that might potentially bind to MTHFD2. We found that there were 30 miRNAs that simultaneously exist in two databases, and were predicted to be implicated in the modulation of MTHFD2 expression level ( Figure 8A). According to the mechanism of miRNA regulating the targeted genes, the predicted miRNA should be down-regulated and negatively correlated with MTHFD2 expression in bladder cancer and miR-383-5p met the two required conditions. Previous studies also showed that miR-383-5p suppressed the proliferation and induced the apoptosis in gastric cancer cells [40]. Hence, miR-383-5p was assumed to be the most potential upstream miRNA of MTHFD2.
Circular RNA (circRNA) is usually considered to be the upstream molecule of miRNA and plays a vital role in regulating miRNA. The data from Starbase database predicted there were the 12 circRNAs that were likely to bind to miRNAs ( Figure 9A). The potential circRNAs involved in miR-383-5p/MTHFD2 axis should be up-regulated and associated with unfavorable prognosis in bladder cancer following the ceRNA hypothesis [41]. The results indicated that, among the most potential up-regulated circRNAs, hsa_circ_0046140 and hsa_circ_0006769, were likely to be implicated in miR-383-5p/MTHFD2 axis ( Table   2). In other words, hsa_circ_0046140 and hsa_circ_0006769/miR-383-5p/MTHFD2 axis were de ned as the possible regulatory pathways in bladder cancer, which is required to be further experimentally veri ed.
The TME is mainly composed of tumor cells, stromal cells and numerous extracellular components (such as cytokines, growth factors, hormones, etc.) [42,43]. It is acknowledged that the immune cell in ltration in TME could inhibit proliferation of tumor cells [44]. For example, cytotoxic CD8 + T cells exert an inhibitory effect on tumor survival by recognizing tumor cell surface antigens, releasing cytotoxins or activating death ligand [45]. We found that the expression of MTHFD2 was positively correlated with the tumor in ltrating immune cells in bladder cancer, by analyzing the immune cell in ltration landscape ( Figure 10A-B).
Immunotherapy, as a supplement to the classical anticancer therapy, has been applied in clinical practice in recent years. It could induce a potent antitumor response through activating the immune system and reverse the immune tolerance, which greatly prolonged the survival of tumor patients and provided a new approach for cancer therapy [46]. Immune checkpoint inhibitors, such as PD-1/PD-L1 and CTLA-4 were widely used in tumor immunotherapy as crucial therapeutic targets. The overexpression of immune checkpoints is an additional factor that determines the e cacy of immunotherapy [47]. We then evaluated the association between MTHFD2 and immune checkpoints and the ndings manifested a positive correlation between MTHFD2 expression and PD-1/PD-L1 and CTLA-4 in bladder cancer. (Figure  11), suggesting that targeting MTHFD2 may help improve the e cacy of immunotherapy for bladder cancer.
However, some recent evidences have found that the immune cell in ltration in TME could prevent tumor cells from being killed. Tumor-associated macrophages (TAMs) have been veri ed to enhance the growth, angiogenesis and invasion of tumor cells by secreting growth factors, cytokines and proteases, and develop resistance to treatment [48]. TAMs also induce chronic in ammation by releasing various cytokines such as tumor necrosis factor, IL-6, IL-12 and transforming growth factor-β, leading to immuneassociated tumor cytotoxicity suppression events [49]. Therefore, the role of activation and in ltration of immune cells in the TME in tumorigenesis, progression and therapeutic effect is still required to be further elucidated.

Conclusion
In summary, MTHFD2 is widely overexpressed in pan-cancer and higher expression of MTHFD2 predicted a worse prognosis, especially in bladder cancer. Hsa_circ_0046140 and hsa_circ_0006769/miR-383-5p/MTHFD2 axis might be involved in multiple biological process concerning tumor proliferation and invasion. Moreover, higher expression of MTHFD2 was correlated with tumor in ltrating immune cells in the TME. Taken together, these ndings demonstrated that targeting MTHFD2 could be a novel approach for bladder cancer therapy.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Data Availability Statement
The data used to support the ndings of this study are included within the manuscript.

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