Integrated bioinformatics analysis identified MTHFD1L as a potential biomarker and correlated with immune infiltrates in hepatocellular carcinoma

Biomarkers are urgently needed for the diagnosis of cancers, such as gastric cancer [1] and liver cancer [2]. Primary liver cancer is one of the most common cause of cancer-related deaths, and is one of the few neoplasms with a steady increasing incidence and mortality in the world [3–5]. There are several main risk factors for liver cancer include excessive alcohol consumption and smoking tobacco, aflatoxin exposure, infection with hepatitis B or C, obesity and diabetes mellitus [3,6–9]. Liver cancer also includes liver hepatocellular carcinoma (LIHC), cholangiocarcinoma and mixed hepatocellular carcinoma, in which LIHC accounts for above 90% of liver cancers [3,10,11]. The prognosis of LIHC depends on the tumor stage at the time of diagnosis, and early detection of liver cancer is easier to obtain treatment options [12].

Therefore, it is urgent to find new and effective biomarkers for LIHC to improve clinical efficacy. Methylenetetrahydrofolate dehydrogenase 1-like (MTHFD1L), is an enzyme involved in tetrahydrofolate (THF) synthesis in mitochondria [13], participates in folic acid cycle and to form formate [14]. Folic acid disorders can contribute to predisposition to some diseases, including immune dysfunction and cancer [15]. Despite the importance of MTHFD1L in the folic acid cycle, knowledge about its role in human diseases, especially cancer, is still largely unknown. A recent study reported that MTHFD1L plays an important role in bladder cancer cell proliferation, colony formation and invasion, which was related to the prognosis [16]. Another study reported that MTHFD1L was involved in colorectal cancer progression, and blocking MTHFD1L reduces the growth of colon cancer cells, thus providing an avenue to target this enzyme with small molecule inhibitors [14]. However, the potential prognostic value of MTHFD1L and its impact on the immune infiltration of LIHC tumors are still unclear.

Plenty of studies have confirmed that the tumor microenvironment (TME) plays a vital role in the occurrence, development and progression of tumors [17–19]. TME refers to the environment surrounding the tumor, including surrounding blood vessels, tumor-infiltrating immune cells (TIICs), fibroblasts, signaling molecules and extracellular matrix [20,21]. TIICs play an important role in the TME and their activation status is different in TME that affects the therapeutic effect of tumor [20]. In the malignant cells inside the tumor, the abundance of TIICs has different degrees of heterogeneity due to the different TMEs, enhancing the immune cell function of tumors is still the main challenge of cancer therapy [22,23].

In this work, we analyzed the RNA-Seq and clinical data of 370 LIHC patients from The Cancer Genome Atlas (TCGA), to identify the relationship between MTHFD1L and clinical information, the expression level of MTHFD1L in liver cancer was verified by Oncomine database and Immunohistochemistry (IHC) study. Furthermore, the data on immune system infiltrates of LIHC from TCGA database were analyzed through the ImmuneCellAI, CIBERSORT and TIMER. The identification of prognostic MTHFD1L suggests the potential role of MTHFD1L in LIHC pathogenesis and progression.

A total of 370 LIHC patients and 51 control samples data were downloaded from the TCGA database (https://portal.gdc.cancer.gov/). Thirty-five tumor and 10 control samples were obtained from Oncomine (http://www.oncomine.org) database [24], a web‐based data mining platform, harbored the microarray databases of most human cancers. Some samples with incomplete clinical information were deleted when LIHC samples were analyzed by comprehensive bioinformatics. A total of 366 patients with RNA sequencing data and clinical information was obtained from TCGA database. The statistical difference of these features between high- and low-expression MTHFD1L groups was evaluated using chi-square test.

According to the level of MTHFD1L expression, we divided 366 LIHC samples into high-expression group and low-expression group, and chi-square test was used to evaluate the correlation between MTHFD1L expression level and clinical characteristics. Then we performed Kaplan–Meier analysis on 366 liver cancer patients through SPSS software (version 24.0, Chicago, IL, United States). A univariate Cox proportional hazard regression analysis was used to evaluate the association between overall survival (OS) time and clinical characteristics in the test cohort. The prognostic value of was considered statistically significant as the P-value was <0.05. Subsequently, a multivariate Cox proportional hazards regression analysis was performed to determine the clinical features than can be used as prognosis.

The composition and abundance of immune cells in the TME have a great influence on tumor progression and the effect of immunotherapy. GraphPad Prism software (version 8.4.2, https://www.graphpad.com/updates/) was used to analyze the correlation between MTHFD1L expression and TIICs, the cutoff was P-value <0.01. The Immune Cell Abundance Identifier (ImmuCellAI) detected the abundance of immune cells, and the differences in the infiltration immune cells of 366 LIHC were analyzed. ImmuneCellAI is a tool for estimating the abundance of 24 immune cells from gene expression datasets, including RNA-Seq and microarray data, of which 24 immune cells are composed of 6 immune cell types and 18 T-cell subtypes (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) [25]. The TIMER (https://cistrome.shinyapps.io/timer/) [26] and CIBERSORT (http://cibersort.stanford.edu/) [27] were used to evaluate the results of tumor immune infiltration analysis in the ImmuneCellAI database.

The protein levels of MTHFD1L (Abcam) in formalin-fixed paraffin embedded liver cancer tissues and adjacent non-tumor tissues were detected by IHC. Briefly, the sections were incubated with primary antibody and secondary antibody, visualized with 3,3′-diaminobenzidine and then counterstained with Hematoxylin for examination by the microscope. Fifteen liver cancer tissues and the paired adjacent normal tissues were chosen randomly. The sections were photographed under a microscope and at least six fields were examined.

All statistical analyses were enforced with corresponding packages of R (Version 3.6.3) and analyzed using SPSS (version 24.0). Kaplan–Meier survival analysis and chi square test was conducted using SPSS (version 24.0). To calculate the 95% CI and HR, we used both the univariate and multivariate models of the Cox analysis. Univariate and multivariate Cox regression analyses were used to evaluate MTHFD1L expression between clinical factors and OS. For all analysis, the threshold P-value <0.05 was considered significant statistically.

In order to implore the MTHFD1L expression levels in LIHC and normal liver tissue samples, the expression level of MTHFD1L in tumors and normal tissue of LIHC were analyzed on TCGA and Oncomine databases. We download the data of LIHC and control samples from the TCGA database, which include 370 LIHC and 51 normal samples. The result showed that MTHFD1L was up-regulated in liver cancer compared with control sample (Figure 1A). Also, a total of 35 tumors and 10 controls from the Oncomine database confirmed the high expression of MTHFD1L, as shown in Figure 1B. To further determine the above results, IHC method was used to detect the expression of MTHFD1L in 15 liver cancer tissues and the paired adjacent normal tissues. We found that MTHFD1L expression was up-regulated in LIHC group (Figure 1C,D).

We next investigated the associations between various clinical parameters and expression of MTHFD1L. The complete clinical information of the LIHC patients included in our study is shown in Table 1. The clinical features include age, gender, vital status, ajcc pathologic m, n and t status, in which we found the low-expression of MTHFD1L group and high-expression of MTHFD1L group were significantly associated with vital status (alive and dead, P-value = 0.038). Next, Kaplan–Meier analysis was explored whether the expression of MTHFD1L has an effect in LIHC OS; the result showed that high expression of MTHFD1L was associated with poor OS (Figure 2B). Additionally, heatmap visually showed the relationship between clinical characteristics with MTHFD1L expression, and the clustering of differentially expressed genes in the high- and low-expression MTHFD1L groups (Figure 2A).

Then, univariate and multivariate Cox regression analyses were performed to research the correlation between OS and clinical characteristics. The univariate Cox regression demonstrated that the ajcc pathologic m, n and t status and level of MTHFD1L were identified as independent risk factors effectively to affect the OS of LIHC patients (Figure 2C). In addition, multivariate Cox proportional hazards regression analysis was performed to determine which clinical features could contribute to the OS. As a result, ajcc pathologic t status and level of MTHFD1L were significantly connected with OS as shown in Figure 2D. These results reveal that the expression of MTHFD1L can be used as an independent clinical factor to affect the OS and prognosis of LIHC patients.

To investigate the effect of MTHFD1L expression on different immune cell types in the LIHC microenvironment, the correlation between MTHFD1L expression and immune score was shown in Figure 3A, the MTHFD1L expression was significantly correlated with immune scores (R = 0.516, P-value <0.001). To gain a comprehensive analysis of immune cells in LIHC, the ImmuCellAI database was utilized to calculate the composition of 24 immune cells in 366 LIHC sample. As was shown in Figure 3B, the immune cells varied significantly among different samples. Therefore, we compared the level of immune cells in the high expression MTHFD1L group and low expression group in the above study. CD8 navie, nTreg, Th17, central memo and neutrophil were negatively correlated with MTHFD1L expression. Cytotoxic, Tr1, Th2, MAIT, dendritic cell (DC), B cell, monocyte, macrophage, NK and CD4 T were positively correlated with MTHFD1L expression. In Figure 3C, immune cells with differences are marked with red squares (P-value <0.01).

Subsequently, we studied the correlation between MTHFD1L expression and the immune cells abundance that are different in Figure 3C based on the Pearson’s correlation coefficient. We found that the high expression of MTHFD1L in Table 2 was positively correlated with the expression of DC, Macrophage and NK cells, but negatively correlated with the expression of CD8 naive cells, and was not significantly correlated with other immune cells (P-value <0.01).

To verify immune cells in LIHC are specifically affected by MTHFD1L expression, we combined three databases CIBERSORT, TIMER and ImmuneCellAI to comprehensively analyze the correlation between various immune cell subtypes and MTHFD1L expression. We found that DCs resting, macrophage M0 and macrophage M2 had a positive correlation with MTHFD1L (Figure 4).

LIHC is one of the most common types of cancer with poor prognosis and high mortality. In recent years, its prevalence has continued to increase [28,29]. Many previous studies have provided information on liver cancer biomarkers [30–32], but so far, the biomarkers that can predict the survival of LIHC patients is still not efficient. Numerous previous studies have reported that MTHFD1L has formyltetrahydrofolate synthase activity and participates in the folic acid cycle, and plays a vital role in supporting tumor growth [14,33,34]. MTHFD1L is overexpressed in multiple solid cancers, such as colorectal cancer [14], tongue squamous cell carcinoma [35], esophageal squamous cell carcinoma [36]. In view of this, we explored the possibility of MTHFD1L as a potential biomarker of LIHC prognosis, and analyzed the interaction between tumor immune infiltration and MTHFD1L expression in LIHC.

Correlation analysis suggested that differences of vital status between high and low expression MTHFD1L groups of LIHC were mainly involved in alive and dead, this illustrates that the expression of MTHFD1L has prognostic value to some extent. An early study identifies MTHFD1L in the folate cycle as an important metabolic pathway in liver cancer cells with the potential for therapeutic [11]. The following Kaplan–Meier analysis also confirmed this view, and the high-expression MTHFD1L group of LIHC had a poor OS curve. Cox regression was performed to identify which clinicopathological features have a significant impact on OS. Considering that the univariate Cox model is insufficient to determine the factors that affect OS, we first used the univariate Cox model to screen for clinical factors related to OS, and the multivariate Cox model to improve the predictive performance of prognostic indicators. Finally, it was found that level of MTHFD1L and ajcc pathologic t status were two risk factors affecting the survival of LIHC patients.

Immune cell infiltration in TME is considered to play an important role in the biological behavior of many kinds of tumors. Until now, many studies have focused on immune landscape and its further clinical application [37–39]. As an essential component of TME, DCs were regarded to be antigen-presenting cells that maintain immune responses. DCs could migrate to tumor tissues through blood circulation in the TME. At different stages of tumor progression, DCs may have different functions, which may be immunostimulatory factors or immunosuppressive factors [40–42]. Macrophages can mediate the proliferation, migration and invasion of tumor cells through a variety of molecular mechanisms, and macrophages infiltration was closely correlated with cancer patient prognosis and therapeutic effect [43,44]. Studies have reported that M0 macrophages can internalize lung tumor-derived exosomes and differentiate into M2 phenotype, while M2 macrophages are characterized by pro-tumorigenic properties [45,46]. Our research provided a comprehensive analysis on the TIICs in LIHC and suggested that M2 macrophages and DCs were associated with the MTHFD1L in LIHC tissue. The analysis results indicated that high expression of MTHFD1L could increase levels of M0 macrophages, M2 macrophages and DCs in the TME of LIHC, which should be helpful for clinical surveillance and treatment of LIHC.

To sum up, MTFHD1L is highly expressed in LIHC through our bioinformatics analyses. Increased MTHFD1L expression correlates with poor prognosis and increased immune infiltration levels in DCs, M0 macrophages and M2 macrophages. Therefore, MTHFD1L may plays an important role in immune cell infiltration and as a prognostic biomarker in LIHC patients.

All the data will be provided on reasonable request from the corresponding author.

The authors declare that there are no competing interests associated with the manuscript.

The authors declare that there are no sources of funding to be acknowledged.

Junhui Chen designed this manuscript. Jie Yang and Qingchun Xu performed experiments and wrote the manuscript. Zhenyu Wang, Jun Wu, Liukui Pan and Kai Huang analyzed the data. Chen Wang edited the manuscript.

The research protocol for the present study was approved by the Ethics Committee of the Second People’s Hospital of Wuhu (number 2018021).

CI

Confidence Interval

DC

dendritic cell

HR

Hazard Ratio

IHC

immunohistochemistry

ImmuCellAI

Immune Cell Abundance Identifier

LIHC

liver hepatocellular carcinoma

MAIT

mucosal-associated invariant T

MTHFD1L

methylenetetrahydrofolate dehydrogenase 1-like

NK

natural killer cell

nTreg

naturally-occurring regulatory T cells

OS

overall survival

RNA-Seq

RNA-Sequence

TCGA

The Cancer Genome Atlas

Th17

T helper cell 17

TIIC

tumor-infiltrating immune cell

TME

tumor microenvironment

Tr1

type 1 T regulatory cells