3.1 Expression of MTHFR across tissues, cancers and pathological stages
We firstly analyze the expression level of gene and protein in cancer to show whether the MTHFR could influence the mechanism of cancer. Expression alterations play important roles in cell function and regulating up and down genes’ transcriptome features. Overall, MTHFR gene expresses significantly different in 24 cancers types. MTHFR gene shows significantly lower expression in 17 tumors, including BRCA, CESC, LUAD, et al. However, the expression level of MTHFR gene has higher expression in 7 tumor samples (Fig. 1A) and exhibited widely expression in normal tissues, such as ALL, GBML, LAML, et al. Specifically, MTHFR shows highest expression in the LAML, lowest expression in the LIHC (Fig. 1A). According to the CPTAC datasets, the total protein level of MTHFR is significantly higher in KIRC and PAAD (Fig. 1B), lower in LUAD. By using “Pathological Stage Plot” module of GEPIA2, the lower MTHFR expression significantly is significantly correlated with the higher pathological stages of cancers in ACC, BLCA, BRCA, et al (Fig. 1C).
3.2 Prognostic significance of MTHFR across cancer types
Tumor samples are divided into high and low expression groups depending on the expression levels of MTHFR. The correlations between MTHFR expression and prognostic significance in different types of cancer are examined. In Fig. 2A, high expressed MTHFR is linked to better prognosis of OS in LGG, LAML, ACC, but associated with the poorer OS prognosis in 7 cancer types, including KIPAN, HNSC, PAAD, et al. According to DSS analysis data (Fig. 2B), high MTHFR expression is predictor of favorable DSS among patients with LGG, UCEC, LUSC, et al., but go against DSS in KIRC, KIRP, GBMLGG, HNSC, PAAD. Additionally, high expression of the MTHFR is related to better PFI in LGG, LUSC, ACC, KICH (Fig. 2C). In contrast, low MTHFR expression could predict the favorable PFI in KIRC, KIRP, HNSC, PAAD. Expression of MTHFR of all cancer types shows no relationship to DFI (Fig. 2D).
Using the Kaplan-Meier plotter tool, analysis of the survival probability of OS presents a significant correlation between low expression of MTHFR and favorable OS prognosis in ACC, LIHC, LAML, LGG, KICH, WT. Meanwhile, high expressed MTHFR could predict the OS among patients with KIRC, HNSC, PAAD, MESO (Fig. 2E).
3.3 Alteration of MTHFR gene across tumor types
Gene alterations of MTHFR are totally different across tumor types (Fig. 3A). MTHFR alteration is most prevalent in UCEC with mutation (4.73%, 25 cases), Amplification (0.76%, 4 cases) and Deep Deletion (0.38%, 2 cases). The amplification type of copy-number alterations (CNA) is the primary type in the ESCA and Uterine Carcinosarcoma cases, which shows an alteration frequency of 2.75% (5 cases) and 1.75% (1 case) respectively. It is worth noting that all Cholangiocarcinoma (2.78%), Pheochromocytoma & Paraganglioma (2.25%), Mesothelioma (1.15%) and Testicular Germ Cell (0.67%) Tumors have copy number deletion of MTHFR (Fig. 3A). Further details are provided in Fig. 3B, including the type, site, and case number of the MTHFR alterations. MTHFR missense mutations are the most prevalent alteration. In 3 cases of UCEC, 1 case of GBM and STAD, R183Q/*/L alteration in the domain is able to induce a missense mutation of the MTHFR gene (Fig. 3B). An alteration to MTHFR protein occurs when R (Arginine) is converted to Q (Glutamine)/L (Leucine) at the 183 sites. Furthermore, we observe the significant differences in MTHFR expression between Neutral, loss and gain in 26 tumors. As shown in Fig. 3C, gain mutation of MTHFR shows highest and lowest expression in the ESCA and BLCA respectively, loss mutation of MTHFR shows highest and lowest expression in the ESCA and LIHC respectively.
3.4 Machine Learning Identifies Stemness Features Associated with MTHFR
Because of related with oncogenic dedifferentiation across tumor types according to previous study [21], machine learning identifies mDNAsi (Fig. 4A) and mRNAsi (Fig. 4B) stemness indices of MTHFR that are Scales ranged from 0 (low) to 1 (high) are analyzed to show the initiation and evolution of cancer. A correlation between cancer stemness index and ICBs response may lead to the discovery of potential therapeutic targets [21]. Using DNA methylation levels, the mDNAsi reflects epigenetic stemness characteristics, while using mRNA expression levels mRNAsi reflects the transcriptomic stemness characteristics. As we can see from Fig. 4C-H, the expression of MTHFR could predict the stemness indices of mDNAsi, mRNAsi, DMPsi, ENHsi, EREG-mDNAsi and EREG-mRNAsi in the most cancers. Specifically, we calculated Pearson correlation in 37 tumors. Results suggests that mDNAsi is significantly negatively correlated with MTHFR in 9 tumors, such as COAD, COADREAD, BRCA, et al. Meanwhile, Correlation of mRNAsi and MTHFR shows negative in 18 tumors. High expression of MTHFR predict low DMPsi significantly among patients with LUAD, COAD, BRCA, KIRP, KIPAN, HNSC, TGCT, BLCA. ENHsi is significantly negatively correlated with MTHFR in COAD, COADREAD, BRCA, et al. EREG-mRNAsi is significantly negatively related with MTHFR in 12 tumors. EREG-mDNAsi is significantly negatively correlated with MTHFR in COAD, COADREAD, BRCA, et al.
3.5 Immune infiltration analysis of MTHFR in cancer
Tumor-infiltrating immune cells, as prominent components of the TME, are closely linked to the tumorigenesis, progression and metastasis. In EPIC algorithm (Fig. 5A), MTHFR is positively correlated with B-cells, cancer-associated fibroblasts (CAFs), CD4 + T-cells, CD8 + T-cells in almost tumors. After a series of analyses by XCELL we obtained 67 types of immune cell infiltration scores of 44 tumor types from a total of 10179 tumor samples (Fig. 5B). Subsequently we calculated the Pearson's correlation coefficient between MTHFR and immune cell infiltration score in each tumor using the corr. test function of R software package psych to determine the significantly related immune infiltration score. Finally, the MTHFR expression showed significantly correlated with immune infiltration in all 44 cancer types. In XCELL algorithm of most cancers, MTHFR is positively correlated with the, CD4 + T-cells, CD8 + T cell, CD4 naive T cell, mast cell and osteoblast and others, negatively correlated with CD8 naive T cell, Th1 and Th2 cell, Mucoid Exopolysaccharide, et al. But the results are quite different when it comes to a specific tumor. These phenomena could be a clue to further research about the tumorigenesis and immunotherapy of MTHFR.
Moreover, CAFs in the stroma of the TME were reported to participate in modulating the function of various tumor-infiltrating immune cells. Herein, using the XCELL, MCPCOUNTER, EPIC and TIDE algorithms, we investigated the possibility of a relationship between the level of infiltration of different immune cells and the expression of the MTHFR gene among TCGA cancer types (Fig. 5C). A statistically positive correlation between MTHFR expression and the estimated infiltration value of CAFs is observed in 8 types of tumors, including BRCA, LIHC, LUSC, et al, but negative correlation in ESCA and LUAD. The scatter plot data of the above tumors processed by one algorithm are presented in Fig. 5D. Except for the negative relationship to ESCA, the MTHFR expression is positively correlated with the CAFs infiltration in the other 7 tumors.
3.6 Enrichment analysis of MTHFR across cancer types
To Study the molecular mechanism by which the MTHFR gene contributes to tumorigenesis, we conducted GO and KEGG Pathway enrichment analyses to identify the targeting of the MTHFR expression-correlated genes and MTHFR-binding proteins. by using the STRING tool, 20 MTHFR-binding proteins are MTRR, KIAA1429, PSMA1, PMSB1, NOS1/2/3 and others (Fig. 6A). The KEGG data of MTHFR-binding proteins suggests that MTHFR is most related with Biosynthesis of amino acids, Pentose phosphate pathway, Carbon metabolism, Parkinson disease, et al (Fig. 6B). According to the GO enrichment analysis of Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) data, MTHFR is most related with Methionine biosynthetic process of BP, Proteasome core complex of CC and Transketolase activity of MF (Fig. 6C).
Utilizing the GEPIA2 tool, we combined all TCGA tumor expression data to identify the top 100 genes associated with MTHFR expression. We combined the genes datasets to perform GO and KEGG enrichment analysis. The KEGG data suggests the most related 11 pathways, including Lysine degradation, the NOTCH signal pathway, Hepatitis B and Human papillomavirus which could be influenced by MTHFR in tumor pathogenesis (Fig. 6D).The GO enrichment analysis further indicates that MTHFR plays an important role in nuclear lumen, nuclear part, regulation of RNA metabolic process, regulation of nucleobase-containing compound metabolic process, nucleoplasm and others (Fig. 6E). In Fig. 6F, the MTHFR is most positively correlated with AFF1 (ALF transcription elongation factor 1, R = 0.60), CPLANE2 (ciliogenesis and planar polarity effector complex subunit 2, R = 0.54), CREBRF (CREB3 regulatory factor, R = 0.54), CYB561D1 (cytochrome b561 family member D1, R = 0.59), WDTC1 (WD and tetratricopeptide repeats 1, R = 0.57). Additionally, the corresponding scatter diagrams indicate positive correlations between MTHFR and the above five genes separately (Fig. 6G).
3.7 Heterogeneity and MTHFR
Our research demonstrates MTHFR expression is a strong predictor of heterogeneity after analyzing the relations with biomarkers MSI, TMB, HRD, MATH, LOH, NEO, Ploidy and Purity. Firstly, MTHFR predicts significantly positive relationship to Gene expression Data of MSI [22] in COAD, COADREAD, LUAD, et al, and negative relationship in BRCA, KIPAN, PRAD, et al (Fig. 7A). We observed MTHFR significantly positively correlated with TMB in COAD, COADREAD, ESCA, ACC and KICH, negatively correlated in HNSC. Accurate detection of HRD [23] is clinical relevance as it is indicative of sensitivity to targeted therapy with poly ADP ribose polymerase inhibitors (PARPi) as well as to DNA damaging reagents [24]. HRD shows significantly positively correlated with MTHFR in ACC, GBMLGG, LGG, but negatively correlated with MTHFR in CESC, BRCA, ESCA, et al (Fig. 7C). MATH is a novel, non-biased, quantitative measure to assess intra-tumor heterogeneity based on next-generation sequencing data. In our study, the MATH shows significantly positively correlated with MTHFR only in GBMLGG. in contrast, MATH is negatively correlated with MTHFR in BRCA, COADREAD, KIPAN, and LUSC (Fig. 7D). Correlations between MTHFR and LOH, NEO, Ploidy and Purity are also significant in specific tumors (Fig. 7E-H).
Different detection methods are used for evaluating the same biological effect, such as DNA for determining MSI status or IHC for detecting MMR protein expression [25]. Nine MMR genes, such as MLH1, PMS2, MSH2, MSH6, EPCAM, MLH3, PMS1, EXO1, have been found in the human MMR system [26]. Low expression of MMR genes meaning less protein will be detected by IHC may be seen as negative or dMMR. In Fig. 7I, MTHFR is significantly positively correlated with most MMR genes in most cancer types. Generally, Mutated MTHFR expressing low level of protein may be correlated with dMMR in IHC. This phenomenon need be verified in future research.
3.8 MTHFR and IRGs across cancer types
Our analyses aim to depict the immunological role of MTHFR which is critical in determining the types of cancers benefit from immunotherapy. Correlation analyses demonstrate that the MTHFR is positively correlated with the most expression of immunoinhibitors, immunostimulators, major histocompatibility complex (MHC), chemokines and their receptors in most tumors (Fig. 8A). MTHFR demonstrates merely positive relations to dozens of IRGs, such as CX3CR1, CCR4, TAP2, HLA-E, TGFBR1, IL6R, TNFSF13, et al. However, heat map displays several cancer types such as CHOL and UCS show almost few relationships between MTHFR and the immune related genes. Based on ICBs, We select the somatic genes CD19, CD274, CD80, CD86 as the most representative genes for analysis. The results demonstrate that MTHFR is significantly positively correlated with CD80/CD86 in BRCA − Basal, COAD, HNSC − HPV (−), KICH, LIHC, LUSC, PRAD and SKCM − Primary, with CD274 in COAD, HNSC, LIHC, SKCM and others, with CD19 in LUSC, BRCA − Basal, COAD and others. (Fig. 8B). To be specific, MTHFR is significantly correlated with CD274 in THYM (P = 2.38e-19), LIHC (P = 1.32e-07), COAD (P = 1.47e-15), with CD86 in COAD (P = 3.9e-09), in LUSC (P = 2.3e-07), in BRCA (P = 1.83e-08), with CD80 in COAD (P = 9.76e-08), with CD19 in LUSC (P = 3.34e-07), et al (Fig. 8C).
3.9 MTHFR act as a potential biomarker of ICBs
Area under curve (AUC) shows how well the test separates the responder and non-responder groups. The larger the area under the ROC curve, the more useful is the measurement to predict treatment response. After calculating data of Ontreatment patients (Fig. 9A-F), MTHFR could act as a biomarker with potential clinical utility for all ICBs therapy (AUC = 0.591, P = 0.028), anti-PD1 therapy (AUC = 0.654, P = 0.003) and anti-ATLA4 therapy (AUC = 0.782, P = 0.002). Precisely, MTHFR could act as blockbuster biomarker for ipilimumab therapy (AUC = 0.9, P < 0.0001) and top-quality cancer biomarker for nivolumab therapy (AUC = 0.712, P = 0.006) to melanoma patients. When coming to the data of pretreatment patients (Fig. 9G-L), MTHFR could act as a top-quality biomarker for pembrolizumab therapy (AUC = 0.704, P < 0.0001) to melanoma patients.
3.10 Verification the expression of MTHFR by IHC
For better clinic application of MTHFR to assist treatment and research, we verified the IHC across cancer types. Figure 10A-B indicates the MTHFR IHC summaries of antibodies HPA076180 and HPA077255 separately in different cancer tissues. Most cancers show the strong and medium positive expression of MTHFR. But HPA076180 shows no positive IHC in 6 cancers, including Glioma, liver cancer, carcinoid, testis cancer, endometrial cancer and skin cancer patients. HPA077255 shows no positive IHC in 5 cancers, including lung cancer, stomach cancer, pancreatic cancer, endometrial cancer and skin cancer. Figure 10C shows the IHC images of positive and negative control of 6 prevalent cancers typically, including colon adenocarcinoma, malignant lymphoma (non-Hodgkin’s type), breast duct carcinoma, cervix squamous cell carcinoma, liver carcinoma hepatocellular and lung adenocarcinoma.