MTHFD2 expression and its association with patients’ survival in pan-cancer
We initially identified 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 findings demonstrated that MTHFD2 was overexpressed in 18 kinds of tumors compared with corresponding normal tissues (Figure 1B-C).
We further examined the associations between the MTHFD2 expression and the patients’ prognostic by Kaplan-Meier curve. The higher expression of MTHFD2 indicated poorer overall survival (OS) (Figure 2) and progress free interval (PFI) (Figure 3) in most cancer types. Taken together, MTHFD2 could be considered as an unfavorable prognostic predictor in pan-cancer.
Upregulation of MTHFD2 and its clinical significance in bladder cancer
Specifically, 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 findings (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 verified 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 findings 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 down-regulated 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 infiltrating immune cells
Some immunological studies were conducted to assess whether expression of MTHFD2 could affect the immune cell infiltration. We found that higher MTHFD2 expression group exhibited a higher level of immune cells infiltration, 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 infiltration 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 infiltrating 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 findings 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 efficacy 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 findings 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.