1. Introduction
Given the complex mechanisms of the development of cancer, it is important to employ the gene analysis to evaluate relevant mechanisms in different kinds of cancer, which contributes to determining new tumor markers for early screening, diagnosis, or evaluating whether they can be regarded as therapeutic targets. The Cancer Genome Atlas (TCGA) is one of the most successful cancer genomics projects to date. The genomic sequence expression methylation and copy-number change data have been generated, analyzed, and provided by TCGA from more than 11,000 patients representing over 30 different cancers.[1] The Gene Expression Analysis Interactive Analysis (GEPIA2) Web server is a valuable and highly cited resource for TCGA-based gene expression analysis of tumor samples and Genotype-Tissue Expression (GTEx) database corresponding to normal tissue samples. Characterized by 198,619 isoforms and 84 cancer subtypes, GEPIA2 has extended the gene expression quantification from the gene level to the transcriptional level, and it can support the analysis of specific cancer subtypes and the comparison between subtypes.[2] Based on the above databases, it is possible to conduct a generalized cancer analysis.
The gene solute carrier family 16 member 1 (SLC16A1, also known as The monocarboxylate transporter isoform 1) is a proton-related monocarbonate transporter that could catalyze the movement of many monocarbonates, such as lactates and phenylpropanates, across serous membranes. Mutations in this gene are correlated with the defect in red blood cell lactate transport. Studies have reported that the mechanism of SLC16A1 exists in different tumors, such as non-small cell lung cancer, bladder cancer, osteosarcoma, breast cancer, and renal carcinoma.[3–8] There are also reports suggesting that SLC16A1 played a key role in promoting the development and progression of tumors by altering cellular metabolic procedures.[9,10]
In this study, TCGA and the Gene Expression Omnibus (GEO) were applied to the pan-cancer analysis of SLC16A1 for the first time. The potential molecular mechanism of SLC16A1 in different tumor genesis and clinical prognosis was also investigated from the aspects of gene expression in different cancers. Besides, we comprehensively analyzed the correlation between gene expression of SLC16A1 and survival prognosis, DNA methylation, gene mutation, protein phosphorylation sites, and immune infiltration analysis, and conducted enrichment analysis of related proteins.
2. Materials and methods
2.1. Gene expression analysis
In the “GENE_DE” module of TIMER2 (Tumor immune estimation resource, version 2) database (TIMER2.0, http://TIMER2.cistrome.org), we studied the differential expression of SLC16A1 between tumor and corresponding normal tissues across all TCGA tumors. For the tumor types with no normal tissues in the database, the “Expression DIY” section of online database GEPIA2 (GEPIA 2, http://GEPIA.cancer-pku.cn) was used to conduct differential expression analysis in conjunction with TCGA and GTEx databases. The “Stage Plot” under this panel can also assess the gene expression levels of tumors in different clinical stages.
2.2. Survival prognosis analysis
The “survival map” region in the “survival analysis” section of GEPIA2 database was utilized to analyze the prognosis of the SLC16A1 gene in 33 tumors. The cutoff values of high (50%) and low (50%) were invoked as the expression thresholds to divide the high and low expression thresholds. Survival bars were clarified with prognostic indicators such as overall survival (OS) and disease-free survival (RFS). Similarly, the survival curve of a single tumor in TCGA can be obtained by selecting different cancer species under this module. To make Kaplan–Meier diagram, we used the Kaplan–Meier plotter (https://kmplot.com/) to analyze the prognosis of lung cancer, gastric cancer, breast cancer, and colon cancer. The prognostic index included OS, RFS, survival without distant metastasis, disease-specific survival (DSS), survival after progression (PPS), survival after the first progression (FP), and progression-free survival (PFS). The above tumor cases were divided into 2 groups by setting the “automatic selection of optimal cutoff point,” and the hazard ratio (HR) and 95% confidence interval were calculated to obtain Kaplan–Meier survival charts.
2.3. Genetic mutation analysis
After logining cBioPortal website (https://www.cbioportal.org/), we chose “TCGA PanCancer Altas Studies” and inputted “SLC16A1” for rapid analysis. In the “Cancer Type Summary” module, different mutation types and mutation frequencies of SLC16A1 gene in 33 tumors were obtained. The mutation site of SLC16A1 gene can be obtained by selecting “mutation,” and the common tumor type on the site can be viewed by clicking on a site. “View 3D Structure” in the lower right corner was used to obtain the 3D Structure of the mutation site information. We commonly used the “Comparision/Survival” panel to obtain the survival curve graphs of mutated and non-mutated states. The prognostic indicators include OS, DSS, disease free survival (DFS), and PFS.
2.4. DNA methylation analysis
We inputted gene SLC16A1 in MEXPRESS website (https://mexpress.be) and selected the Brain Lower Grade Glioma (LGG) in TCGA database to analyze the DNA methylation level of SLC16A1 of multiple probes (CG07176692, CG19645639, CG07967156, etc.) The sample β value, Benjamini-Hochberg adjusted P value and Pearson correlation coefficient (R) value were also obtained. LGG of nervous system tumors was selected and the SLC16A1 was inputted in DiseaseMeth2.0 (the 2.0 version of human disease methylation database). We selected “450K based on array technology” and chose gene sequences NM_003051 and NM_001166496. To carry out t test, we acquired the difference of methylation expression between normal tissue and tumor tissue. Methylation of SLC16A1 in LGG cancer species was analyzed using UALCAN website, and the difference in methylation expression between mutation and non-mutation was obtained.
2.5. Protein phosphorylation analysis
We also used UALCAN cancer database (ualcan.path.uab.edu/home) in the “CTPAC analysis” to analyze the protein phosphorylation of SLC16A1. Total protein expression levels were analyzed in normal tissues and renal clear cell carcinoma, lung cancer, endometrial carcinoma, breast cancer, ovarian cancer, and colon cancer. Differences of the expression level of phosphorylation sites (including Np_001159968. 1: S467, S498, S461, S483, T466S467) between normal tissue and tumor tissue in breast cancer, ovarian cancer, renal clear cell cancer, endometrial cancer, and pediatric brain cancer can be obtained by “CTPAC Analysis” from the UALCAN database.
2.6. Immune infiltration analysis
We applied the “immune module” of TIMER2 database, inputted SLC16A1, and selected CD8+ T-cells, cancer-associated fibroblasts in the immune infiltration factor. The relationships between SLC16A1 gene expression and CD8+ T cells as well as fibroblasts were obtained by CIBERSORT, Cibersort-ABS, QUANTISEQ, XCELL, and other algorithms. The data could be visualized by clicking different color blocks, and the scatter plot could be obtained. In addition, we analyzed the correlation of SLC16A1 expression with tumor mutation burden (TMB) and microsatellite instability (MSI) across all tumors of TCGA by Spearman’s correlation.
2.7. SLC16A1-related gene enrichment analysis
First, the SLC16A1 binding protein of “Homo Sapiens” was screened out using the STRING database (string-db.org), and the conditions were set as followed: The network type = “full network”, the active interaction source = “all options,” the meaning of the network edge = “evidence,” the minimum interaction score = “low confidence (0.150),” and the maximum number of interactors displayed = “no more than 50 interactors in 1st shell.” We then used “Similar Genes Detection” of GEPIA2 to obtain the top 100 SLC16A1-related target Genes. We randomly selected several related genes to conduct correlation analysis under the “Correlation Analysis” module and obtained the correlation scatter plot. The point plot was log2 TPM, indicating the P value and the correlation coefficient (r). These genes were input into the exploration module of the TIMER2 database, and the gene correlation heat map was obtained. In addition, we combined these two sets of data for KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis. We uploaded this list of 151 genes on the DAVID (Database for annotation, visualization, and integrated discovery) (ncifcrf.gov). The data with P value less than .05 were screened for the KEGG bubble map using “R studio” software (version 4.0.3, https://www.r-project.org/), and the pathway function annotation with statistical significance was obtained. The GO_MF (molecular function of gene body) enrichment analysis was executed for 151 genes, and the bubble map was obtained.
2.8. Statistical analysis
R software version 4.0.0 and SPSS version 24.0 were used for statistical analyses. Two-sided P value < .05 was considered statistically significant.
3. Results
3.1. Differential expression of SLC16A1 in tumors and normal tissues
In our research, we aimed at exploring the oncogenic role of SLC16A1. TIMER2 was utilized to analyze the expression of SLC16A1 in different cancers in TCGA database. Compared with the normal control group, the expression of SLC16A1 was markedly increased in esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), renal chromatic cell carcinoma (KICH), renal clear cell carcinoma (KIRC), lung squamous cell carcinoma (LUSC), prostate cancer, and gastric cancer (STAD). The expression levels of SLC16A1 in Breast invasive carcinoma (BRCA), Cholangiocarcinoma (CHOL), Colon adenocarcinoma (COAD), and Uterine Corpus Endometrial Carcinoma were lower than those in normal tissues (Fig. 1A). After adding normal tissues from the GTEx dataset, a further study was performed on the differences in SLC16A1 expression between normal and tumor tissues of CESC (Cervical squamous cell carcinoma and endocervical adenocarcinoma), DLBC (Lymphoid Neoplasm Diffuse Large B-cell Lymphoma), ESCA, GBM, HNSC, KICH, KIRC, LGG, LUSC, PAAD (Pancreatic adenocarcinoma), SKCM (Skin Cutaneous Melanoma), STAD, and THYM (Thymoma). Compared with the normal control group, SLC16A1 was highly expressed in these cancer species with a statistically significant difference (Fig. 1B). In addition, the results of CPTAC dataset showed that the expression of SLC16A1 protein was higher in tumor tissues of renal clear cell carcinoma, lung cancer, and endometrial carcinoma than that of normal tissues, and was lower in breast cancer, ovarian cancer, and colon cancer (Fig. 1C) (Table S1, Supplemental Digital Content, https://links.lww.com/MD/I652). Furthermore, the “pathological staging map” module of HEPIA2 was employed to explore the correlation between the expression of SLC16A1 and the pathological staging of cancer. The results demonstrated that there was a statistical correlation in Bladder Urothelial Carcinoma (BLCA), BRCA, Kidney renal papillary cell carcinoma (KIRP), testicular germ cell tumors (TGCT), KICH, and PAAD (Fig. 1D).
;)
Figure 1.: Expression level of SLC16A1 gene in different tumor and pathological stages. (A) The expression level of SLC16A1 gene in TCGA tumors and corresponding normal tissues through TIMER2. (B) Difference of SLC16A1 expression between DLBC, LAML, LGG, SKCM, THYM, and the corresponding normal tissues of GTEx project. (C) The difference of SLC16A1 total protein expression levels in UCEC, COAD, and RCC. (D)The expression levels of SLC16A1 in different pathological stages (stage I, stage II, stage III and stage IV) analyzed in TCGA database. ***P < .001, **P < .01, *P < .05. BLCA = Bladder Urothelial Carcinoma, BRCA = breast invasive carcinoma, COAD = colon adenocarcinoma, DLBC = Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, GTEx = Genotype-Tissue Expression, KICH = renal chromatic cell carcinoma, LAML = Acute Myeloid Leukemia, LGG = Brain Lower Grade Glioma, RCC = renal clear cell carcinoma, SKCM = skin cutaneous melanoma, SLC16A1 = the solute carrier family 16 member 1, TCGA = The Cancer Genome Atlas, THYM = thymoma.
3.2. The expression level of SLC16A1 is correlated with the prognosis of patients with different tumors
According to the expression level of SLC16A1, the cases with tumors were divided into the high expression group and the low expression group. GEPIA2 database and Kaplan–Meier Plotter website were employed to explore the relationship between the expression level of SLC16A1 and the prognosis of different tumors. GEPIA2 database analysis of the pan-cancer cohorts showed that the high expression of SCL16A1 was correlated with the poor OS of ACC (P = .025), BLCA (P = .012), KICH (P = .001), and LGG (P < .01), PAAD (P = .032), pheochromocytoma and paraganglioma (PCPG) (P = .028), and SKCM (P = .013) (Fig. 2A). As the results of DFS analysis, the high expression of SLC16A1 was correlated with the poor prognosis in TCGA cases with ACC (P = .004), KICH (P = .019), LGG (P < .01), ovarian serous cystadenocarcinoma (OV) (P = .018) and STAD (P = .013) (Fig. 2B). Furthermore, integrated analysis of GEO, EGA, and TCGA databases in the Kaplan–Meier Plotter website reveled that the high expression of SLC16A1 was correlated with the poor RFS (P = .002), OS (P = .012), and PPS (P = .046) in breast cancer cases. Similarly, in ovarian cancer, the high expression of SLC16A1 was correlated with poor RFS (P < .01) and OS (P = .013). Moreover, the high expression of SLC16A1 was correlated with poor OS (P < .01) and FP (P < .01) in lung cancer. In contrast, the high expression of SLC16A1 was associated with favorable prognostic OS (P < .01), FP (P = .002), and PPS (P < .01) in gastric cancer. There was no statistically significant association between the expression of SLC16A1 and the OS, PFS, RFS, and DSS in liver cancer. The expression of SLC16A1 was correlated with the prognosis of breast cancer, ovarian cancer, lung cancer, and stomach cancer, but was not correlated with the prognosis of liver cancer (Fig. 2C). The above results demonstrated that there were differences between the expression of SLC16A1 and the prognosis in different tumors.
Figure 2.: Relationship between SLC16A1 gene expression and survival prognosis of TCGA tumor. GEPIA2 was applied to detect the correlation of SLC16A1 gene expression with OS (A) and DFS (B) in different tumors. (C) Survival maps were drawn using Kaplan–Meier plotter to study the correlation between gene expression and survival prognosis of lung cancer, ovarian cancer, lung cancer, stomach cancer and liver cancer. BLCA = Bladder Urothelial Carcinoma, BRCA = breast invasive carcinoma, CESC = cervical squamous cell carcinoma and endocervical adenocarcinoma, CHOL = cholangiocarcinoma, COAD = colon adenocarcinoma, DFS = disease free survival, DLBC = Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, ESCA = esophageal carcinoma, GBM = glioblastoma multiforme, GEPIA2 = Gene Expression Analysis Interactive Analysis, HNSC = head and neck squamous cell carcinoma, KICH = renal chromatic cell carcinoma, KIRC = renal clear cell carcinoma, KIRP = kidney renal papillary cell carcinoma, LGG = Brain Lower Grade Glioma, LUSC = lung squamous cell carcinoma, OS = overall survival, OV = ovarian serous cystadenocarcinoma, PAAD = pancreatic adenocarcinoma, PCPG = pheochromocytoma and paraganglioma, PRAD = prostate cancer, SARC = sarcoma, SKCM = skin cutaneous melanoma, SLC16A1 = the solute carrier family 16 member 1, STAD = gastric cancer, TGCT = testicular germ cell tumors, THCA = thyroid carcinoma, THYM = thymoma, UCEC = Uterine Corpus Endometrial Carcinoma.
3.3. Correlation between the DNA methylation level and the gene expression of probes at promoter region
MEXPRESS was adopted to explore the DNA promoter methylation level of SLC16A1 in LGG cancer. It was found that the gene expression of the probe of promoter regions showed a significant negative correlation with methylation. For example, CG07176692 (r = −0.181, P = .001) and CG19645639 (r = −0.156, P = .001), instead of the probe of non-promoter regions such as CG18130079 (r = 0.126, P = .01), CG24738650 (r = 0.095, P = .05), CG01792145 (r = 0.157, P = 0. 001), CG07572909 (r = 0.266, P = .001), CG07967176 (r = 0.309, P = .001), and CG20614262 (r = 0. 171, P = .001) showed a significant positive correlation with methylation (Fig. 3A). Based on the analysis of DiseaseMeth2.0, the methylation of SLC16A1 in normal brain tissues was lower than that of low-grade glioma tissues (P < .01) (Fig. 3B) (Table S2, Supplemental Digital Content, https://links.lww.com/MD/I653). According to UALCAN analysis, the level of methylation in LGG tissues with TP53 mutation was higher than that of no methylation (P < .01) (Fig. 3C).
Figure 3.: DNA methylation analysis of SLC16A1 gene. (A) MEXPRESS method was used to analyze the DNA methylation levels of multiple probes in SLC16A1 in low-grade gliomas. (B) The difference of methylation expression between normal tissue and LGG tumor tissue was compared. (C) Promoter methylation level of SLC16A1 with or without TP53 mutation in LGG. LGG = Brain Lower Grade Glioma, SLC16A1 = the solute carrier family 16 member 1.
3.4. Analysis of SLC16A1 gene mutation among different cancers
By analyzing the genetic mutation status of SLC16A1 in different cancer species in TCGA database, it was found that the “mutations” (>3.5%) type was the most frequent mutation of SLC16A1, which were also the main type of mutation in endometrial cancer, skin melanoma, and lung squamous cell cancer. “Depth loss” is the main mutation in pheochromocytoma and paraganglioma cases, with a change frequency of about 3.4% (Fig. 4A). The type, site, and number of cases with SLC16A1 mutation were further derived through the cBioPortal tool. The 3D structure of the mutation site information was also obtained (Fig. 4B). It was found that the “missense” mutation of SLC16A1 was the main type of genetic alteration. The mutation site R328Q/* in the MFS_1 domain of SLC16A1 gene was the site with the highest mutation frequency. R328Q/* alteration in the MFS_1 domain, which was detected in 3 cases with uterine endometrioid carcinoma, 1 case with colon adenocarcinoma, 1 case with rectal adenocarcinoma and 1 case with HNSC, is able to induce a missense mutation of the SLC16A1 gene, and a nonsense mutation of this gene was also found in 1 case with stomach adenocarcinoma (Fig. 4C). In addition, we further investigated the potential association between SLC16A1 mutation and clinical survival in patients with different types of cancer. The results demonstrated that there was no statistically significant improvement in the prognosis (OS, DSS, DFS, and PFS) of those cases with SLC16A1 mutation, compared with those without SLC16A1 mutation (Fig. 4D).
Figure 4.: The type of SLC16A1 gene mutation in TCGA tumors was correlated with prognosis. (A) Different types and frequencies of gene mutations in TCGA tumors. (B) The 3D structure showed the protein domain and mutation site of SLC16A1 gene. (C) The mutation sites with the highest mutation frequency in SLC16A1 gene and the tumor cases where mutations at this site can be detected. (D) The difference of survival and prognosis of tumor with or without gene mutation. COAD = colon adenocarcinoma, HNSC = head and neck squamous cell carcinoma, SLC16A1 = the solute carrier family 16 member 1, STAD = gastric cancer, TCGA = The Cancer Genome Atlas, UCEC = Uterine Corpus Endometrial Carcinoma.
3.5. The difference in SLC16A1 phosphorylation level between normal and primary tumor tissues
Considering the importance of protein phosphorylation in post-translational modification, we checked the phosphorylation level of SLC16A1 between normal and primary tumor tissues. The phosphorylation level at S467 of SLC16A1 protein was significantly lower in normal tissues than that of ovarian cancer (P < .01) (Fig. 5B) and renal clear cell carcinoma (P < .01) (Fig. 5D). In colon cancer (P < .01) (Fig. 5C) and breast cancer (P < .01) (Fig. 5A), the phosphorylation level of S467 decreased. The expression of T478 phosphorylation site in colon cancer was lower than that in normal tissues (P < .01) (Fig. 5C). The expression of S498 phosphorylation site in ovarian cancer was lower than that in normal tissues (P < .01) (Fig. 5B). It is necessary to perform further analysis to explore the potential role of phosphorylation sites of SLC16A1 in tumorigenesis.
Figure 5.: Analysis of SLC16A1 phosphorylation site expression in different tumors and normal tissues. (A) Breast cancer. (B) Ovarian cancer. (C) Colon cancer. (D) RCC. (E) UCEC. RCC = renal clear cell carcinoma, SLC16A1 = the solute carrier family 16 member 1, UCEC = Uterine Corpus Endometrial Carcinoma.
3.6. Correlation between tumor infiltrating immune cells and SLC16A1 expression in different cancer species
As an important part of the tumor microenvironment, tumor-infiltrating immune cells are closely connected with the occurrence, development, and metastasis of cancer. It has been reported that tumor-associated fibroblasts in the tumor microenvironment stroma are involved in regulating the function of various tumor-infiltrating immune cells. Based on most algorithms, a statistically negative correlation could be observed between CD8+ T cell immune infiltration and SLC16A1 expression in HNSC (especially HPV + HNSC), LUSC, SARC (Sarcoma), TGCT. Moreover, we observed a statistical positive correlation in the tumors of THYM (Fig. 6A and B). In addition, it could be observed that the SLC16A1 mRNA was positively correlated with the invasion value of cancer-related fibroblast in BLCA, BRCA, KIRC, KIRP, PAAD, PCPG, and thyroid carcinoma (Fig. 6C and D).
Figure 6.: Immune infiltration analysis. (A) Association between SLC16A1 gene expression and immune infiltration of CD8+ T cells, (B) Visualized the correlation and obtained the scatter diagram. (C) Correlation of SLC16A1 gene expression and tumor-associated fibroblasts in TCGA tumors. (D) Visualization and obtained the scatter diagram. ACC = adrenocortical carcinoma, BLCA = bladder urothelial carcinoma, BRCA = breast invasive carcinoma, CESC = cervical squamous cell carcinoma and endocervical adenocarcinoma, CHOL = cholangiocarcinoma, COAD = colon adenocarcinoma, DLBC = lymphoid neoplasm diffuse large B-cell lymphoma, ESCA = esophageal carcinoma, GBM = glioblastoma multiforme, HNSC = head and neck squamous cell carcinoma, KICH = renal chromatic cell carcinoma, KIRC = renal clear cell carcinoma, KIRP = kidney renal papillary cell carcinoma, LGG = brain lower grade glioma, LIHC = liver hepatocellular carcinoma, LUAD = lung adenocarcinoma, LUSC = lung squamous cell carcinoma, MESO = mesothelioma, OV = ovarian serous cystadenocarcinoma, PAAD = pancreatic adenocarcinoma, PCPG = pheochromocytoma and paraganglioma, PRAD = prostate cancer, READ = rectum adenocarcinoma, SARC = sarcoma, SKCM = skin cutaneous melanoma, SLC16A1 = the solute carrier family 16 member 1, STAD = gastric cancer, TCGA = The Cancer Genome Atlas, TGCT = testicular germ cell tumors, THCA = thyroid carcinoma, THYM = thymoma, UCEC = uterine corpus endometrial carcinoma, UCS = uterine carcinosarcoma, UVM = uveal melanoma.
3.7. Correlation analysis of SLC16A1 with TMB and MSI
We studied the correlation of SLC16A1 expression with TMB and MSI in different TCGA tumors. The results indicated that the aberrant expression of SLC16A1 was associated with MSI in the patients with OV (P < .01), UCS (P < .01), KICH (P = .043), DLBC (P < .01), STAD (P < .01), COAD (P = .015), BRCA (P = .004), ESCA (P = .016), and TGCT (P = .006) (Fig. 7A). Similarly, the aberrant expression of SLC16A1 was associated with TMB in the patients with KIRC (P = .002), STAD (P < .01), LIHC (P = .008), KIRP (P = .041), TGCT (P = .006), and LUAD (P = .001) (Fig. 7B).
Figure 7.: The potential correlation between SLC16A1 gene expression and MSI and TMB in TCGA tumors. (A) Radar chart of the correlation of SLC16A1 expression with MSI. (B) Radar chart of the correlation of SLC16A1 expression with TMB. BLCA = Bladder Urothelial Carcinoma, BRCA = breast invasive carcinoma, CESC = cervical squamous cell carcinoma and endocervical adenocarcinoma, CHOL = cholangiocarcinoma, COAD = colon adenocarcinoma, DLBC = Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, ESCA = esophageal carcinoma, GBM = glioblastoma multiforme, HNSC = head and neck squamous cell carcinoma, KICH = renal chromatic cell carcinoma, KIRC = renal clear cell carcinoma, KIRP = kidney renal papillary cell carcinoma, LGG = Brain Lower Grade Glioma, LUSC = lung squamous cell carcinoma, MSI = microsatellite instability, OV = ovarian serous cystadenocarcinoma, PAAD = pancreatic adenocarcinoma, PCPG = pheochromocytoma and paraganglioma, PRAD = prostate cancer, SARC = sarcoma, SKCM = skin cutaneous melanoma, SLC16A1 = the solute carrier family 16 member 1, STAD = gastric cancer, TCGA = The Cancer Genome Atlas, TGCT = testicular germ cell tumors, THCA = thyroid carcinoma, THYM = thymoma, TMB = tumor mutation burden, UCEC = Uterine Corpus Endometrial Carcinoma.
3.8. Enrichment analysis of SLC16A1 related proteins
To further explore the molecular mechanism of SLC16A1 gene in tumorigenesis, SLC16A1 binding-protein and SLC16A1 expression-related genes were screened out through a series of pathway enrichment analyses. A total of 51 SCL16A1 binding proteins were obtained with the STRING tool (Fig. 8A). Besides, the first 100 genes correlated with SCL16A1 expression were obtained through GEPIA2 in combination with all tumor expression data from TCGA. The cross-over analysis showed that the first 100 genes and the top 50 protein-interacting genes had an intersection gene, namely ZNF598 (Fig. 8B). As shown in Figure 8D, the expression level of SCL16A1 was positively correlated with the expression levels of CORO1C (R = 0.44), GNAI3 (R = 0.43), MAP4K4 (R = 0.37), PGM2(R = 0.36), PUS7 (R = 0.35), WDR3 (R = 0.35). The corresponding heat map data also showed that SLC16A1 was positively correlated with the six genes in most specific cancer types (Fig. 8C). In this study, the two data sets were combined for KEGG and GO enrichment analysis. KEGG analysis indicated that the role of SLC16A1 in tumor pathogenesis may be mainly related to “central carbon metabolism in cancer,” “glycolysis and glycometabolic synthesis” and “biosynthesis of antibiotics” (Fig. 8E). GO enrichment analysis further demonstrated that most of the genes were correlated with anion transmembrane transporter activity or secondary active transmembrane transporter activity, such as “glucose transmembrane transporter activity,” “organic anion transmembrane transporter activity,” “cadherin binding and symporter activity,” “monocarboxylic acid transmembrane transporter activity,” and “solute: sodium symporter activity” (Fig. 8F and G).
Figure 8.: Enrichment analysis of SLC16A1 related genes was performed. (A) The protein interaction map of SLC16A1 binding protein. (B) The SLC16A1 binding protein and related genes were cross-analyzed. (C) The first 100 genes related to SLC16A1 were identified in the TCGA database, and the correlation between SLC16A1 gene expression and CORO1C, GNAI3, MAP4K4, PGM2, PUS7, and WDR3 was analyzed. (D) Visualization of correlation analysis. (E) The KEGG pathway of SLC16A1 binding protein and related genes was analyzed, and the bubble map was made. (F) GO analysis of the above genes was performed to make a bubble map. (G) The top five items of the GO analysis bubble chart. BLCA = Bladder Urothelial Carcinoma, BRCA = breast invasive carcinoma, CESC = cervical squamous cell carcinoma and endocervical adenocarcinoma, CHOL = cholangiocarcinoma, COAD = colon adenocarcinoma, DLBC = Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, ESCA = esophageal carcinoma, GBM = glioblastoma multiforme, GO = Gene Ontology, HNSC = head and neck squamous cell carcinoma, KEGG = Kyoto Encyclopedia of Genes and Genomes, KICH = renal chromatic cell carcinoma, KIRC = renal clear cell carcinoma, KIRP = kidney renal papillary cell carcinoma, LGG = Brain Lower Grade Glioma, LUSC = lung squamous cell carcinoma, OV = ovarian serous cystadenocarcinoma, PAAD = pancreatic adenocarcinoma, PCPG = pheochromocytoma and paraganglioma, PRAD = prostate cancer, SARC = sarcoma, SKCM = skin cutaneous melanoma, SLC16A1 = the solute carrier family 16 member 1, STAD = gastric cancer, TCGA = The Cancer Genome Atlas, TGCT = testicular germ cell tumors, THCA = thyroid carcinoma, THYM = thymoma, UCEC = Uterine Corpus Endometrial Carcinoma.
4. Discussion
It has been suggested in some studies that SLC16A1 is highly expressed in most tumors and it is correlated with survival prognosis. However, correlation analysis in different tumors showed different conclusions.[11–13] Although emerging research has evaluated the expression level and function of SLC16A1 in certain tumors, its roles across various tumor types, especially the aspects related to prognostic potential and clinical significance, were not systematically studied. Therefore, SLC16A1 in a total of 33 different tumors was comprehensively detected based on TCGA, CPTAC, and GEO database, as well as the molecular characteristics of gene expression, gene mutation, DNA methylation, and protein phosphorylation. In this study, RNA-seq analysis was performed on the tissue samples from human individuals representing 27 different tissues in order to determine tissue-specificity of all protein-coding genes. It was found that SLC16A1 has the highest expression in the heart, followed by the colon and testis. In fact, SLC16A1 can be expressed in all detected tissues (with the low expression in salivary glands and pancreas), showing a low degree of tissue specificity. We found that SLC16A1 was highly expressed in most tumors, such as ESCA, GBM, HNSC, KICH, KIRC, LUSC, prostate cancer, and STAD. While the expression levels of SLC16A1 were significantly lower in BRCA, CHOL, COAD, and Uterine Corpus Endometrial Carcinoma than those in normal tissues. Therefore, SLC16A1 gene may play different roles in different types of cancer. We also observed a statistical correlation between SLC16A1 expression and the pathological stages in BLCA, KIRP, TGCT, KICH, and PAAD. These results suggested that SLC16A1 may play an oncogenic role in promoting cell proliferation.
It was reported that in previous studies that the expression level of SLC16A1 in tumors is significantly higher than that in normal tissues, which is correlated with poor clinical outcomes in patients with cutaneous malignant melanoma. Instantaneous knockout of SLC16A1 can markedly inhibit the growth and invasion of malignant melanoma cells.[9] Our research confirmed this tendency and found that the high expression of SLC6A1 was correlated with a poor overall survival prognosis of cutaneous melanoma (P = .013). To be noted, the expression level of this gene in gastric cancer was higher than that in normal tissues. However, survival analyses in GEPIA2 and Kaplan–Meier plotter showed opposite results. According to GEPIA2 database, SLC16A1 was correlated with worse DFS, while the Kaplan–Meier plotter survival analysis indicated the high expression of SLC16A1 was associated with longer OS, FP, and PPS. This discrepancy may be attributed to different samples from separate databases. Therefore, more clinical data were required to support the correlation analysis between gene expression and tumor prognosis in the generalized cancer. Besides, further subgroup analysis of different prognostic indexes was required, which would conduce to approaching the clinical studies in the real world.
As for lung cancer, it was found in our study that SLC16A1 was highly expressed in LUAD and LUSC based on TCGA and GEO databases, and the statistically significant prognostic correlation between the high expression of SLC16A1 and the poor prognosis in lung cancer was found by the Kaplan–Meier Plotter. Eilertsen et al[14] proposed that the expression of SLC16A1 in tumor cells and mesenchymal cells had a significant and independent effect on tumor prognosis, but the effect was opposite in these cells. SLC16A1 is an independent positive prognostic factor in cancer cells. In mesenchymal cells, SLC16A1 is an independent negative predictor, and it has been demonstrated to be correlated to different functional mechanisms of genes. Therefore, benign and malignant tumors may have different relevance between gene expression and prognostic even in the same organ.
DNA methylation change has been observed in various cancers and was deemed to a carcinogenic factor. For LGG patients of TCGA, we observed a potential correlation between the reduced DNA methylation status at the promoter region and high expression of SLC16A1. Additional evidence is warranted for the potential tumor-promoting role of SLC16A1 DNA methylation in LGG tumorigenesis.
In this study, CTPAC database was utilized to analyze the phosphorylation sites of SLC16A1 in tumors. It was found that the common sites were S467 and S498. The results also showed that the total protein of SLC16A1 and the phosphorylated SLC16A1 at the site of S467 expressed higher in renal clear cell carcinoma and endometrial carcinoma than in normal tissues, and expressed lower in breast and colon cancer. At present, there is no relevant study elucidating the molecular mechanism of abnormal phosphorylation sites in the process of tumor genesis, proliferation, and metastasis.
The tumor microenvironment is essential for tumorigenesis and tumor development. Hence, exploration was conducted to investigate the relationship between SLC16A1 and the immune infiltrating cells, such as fibroblasts and CT8+ T cells in the tumor microenvironment. Cancer-related fibroblasts is an important component of tumor microenvironment, and its role in tumor development cannot be ignored.[15] Cancer-related fibroblast (CAF) could not only release a large number of cytokines to regulate the biological behavior of tumors, but also mediate the proliferation, invasion, and metastasis of tumor cells by promoting the generation of tumor vessels, the remodeling of tumor stroma, the generation of drug resistance, and immunosuppression.[16] In this study, it can be found that there was a statistically positive correlation between SLC16A1 and CAF in various tumors. These results may indicate that the high expression of SLC16A1 in these tumors, along with the increase of CAF, would further promote the proliferation and metastasis of tumor cells. CD8+ T cells, also known as cytotoxic T cells, can direct contact with target tumor cells and release cytotoxic granules and death ligands to kill the target tumor cells. Meanwhile, these cells may induce the apoptosis of target cells by releasing some cytokines. Memory T cells can also play an important role in immune surveillance against tumor cells. SLC16A1 is negatively correlated with CD8+ T cells in some cancers, which would inhibit its killing effect on tumor cells and lead to a poor prognosis. In tumors with highly expressed SLC16A1 genes, the immune infiltrating cells can inhibit the killing effect on tumor cells, and promote their proliferation and metastasis, resulting in the poor prognosis of patients. MSI and TMB status were prognostic indicators in multiple cancers and were commonly regarded as the notable biomarkers associated with the immunological therapy effect. In this research, we analyzed the correlation between SLC16A1 expression, MSI and TMB in all cancers of TCGA database. Based on the analysis, SLC16A1 had a statistically significant correlation with MSI and TMB in STAD and TGCT. Besides, according to the analysis data of gene expression, SLC16A1 was highly expressed in STAD. Therefore, the perfect immunological examination would contribute to the diagnosis and treatment of STAD, indicating a favorable prognosis of immunotherapy. It has been reported that the monocarboxylic acid transporter SLC16A1 is involved in a series of cell biological processes, such as the absorption of short-chain fatty acids in intestinal microflora,[17] the lactate transport in astrocytes and retinal pigment epithelial cells,[9,18] and the insulin secretion regulation in pancreatic β cells.[19] There was a functional correlation between SLC16A1 and clinical diseases, especially tumors, such as lung cancer, pancreatic cancer, colorectal cancer, oral squamous cell carcinoma, glioblastoma, melanoma.[20–24] Although it is required to be demonstrated whether the pathogenesis of SLC16A1 is the same in different tumors, it was found that the mechanism of the action of this gene in the tumorigenesis is inseparable from glycolysis and glucose metabolism synthesis, which is consistent with the enrichment analysis results in this study.[25] A study by Gray et al revealed that SLC16A1 could regulate the tumor migration by activating the HGF/C-MET pathway, in addition to its function as a proton transporter. HGF/C-MET pathway could induce epithelial-mesenchymal transformation in tumor cells, characterized by the absence of cell-to-cell adhesion, which would cause more active movement, invasion, and metastasis of tumor cells.[26,27] More importantly, over 80 years after Warburg’s groundbreaking research, metabolic targeting remains a largely untapped domain in cancer therapies. Based on the fact that SLC16A1 is a key molecular effector of tumor metabolism, it can be regarded as promising therapeutic target.
However, there are also limitations associated with the present study. Firstly, as mentioned above, the comprehensive pan-cancer analysis was performed using public databases. The data from different databases may derive from different workflow, analysis method and population. These may have inevitably caused bias. This bias may be negligible in further studies with large sample sizes. Secondly, our study mainly determined through the bioinformatic approach. The potential function and mechanisms of SLC16A1 should to be further verified in vivo/in vitro experiments.
5. Conclusions
In summary, our studies, for the first time, showed that the SLC16A1 expression was statistically correlated with survival prognosis, DNA methylation, gene variation, protein phosphorylation, immune infiltration analysis, as well as TMB and MSI. We also carried out an enrichment analysis of SLC16A1 in different cancers. The current study emphasized the critical role of SLC16A1 in diagnosis and prognosis of tumors, and its potential as a promising therapeutic target for tumors.
Author contributions
Data curation: Yang Li.
Funding acquisition: Xinna Deng.
Writing – original draft: Lingyun Chen.
Writing – review & editing: Xinna Deng.
References
[1]. Wang Z, Jensen MA, Zenklusen JC. A Practical guide to the cancer genome atlas (TCGA). Methods Mol Biol. 2016;1418:111–41.
[2]. Tang Z, Kang B, Li C, et al. GEPIA2: an enhanced web server for large-scale expression profiling and interactive analysis. Nucleic Acids Res. 2019;47:W556–60.
[3]. Dai J, Li Z, Amos CI, et al. Systematic analyses of regulatory variants in DNase I hypersensitive sites identified two novel lung cancer susceptibility loci. Carcinogenesis. 2019;40:432–40.
[4]. Afonso J, Santos LL, Miranda-Goncalves V, et al. CD147 and
SLC16A1-potential partners in bladder cancer aggressiveness and cisplatin resistance. Mol Carcinog. 2015;54:1451–66.
[5]. Zhao Z, Wu MS, Zou C, et al. Downregulation of
SLC16A1 inhibits tumor growth, metastasis and enhances chemotherapeutic efficacy in osteosarcoma through regulation of the NF-kappaB pathway. Cancer Lett. 2014;342:150–8.
[6]. Hou L, Zhao Y, Song GQ, et al. Interfering cellular lactate homeostasis overcomes Taxol resistance of breast cancer cells through the microRNA-124-mediated lactate transporter (
SLC16A1) inhibition. Cancer Cell Int. 2019;19:193.
[7]. Guo C, Huang T, Wang QH, et al. Monocarboxylate transporter 1 and monocarboxylate transporter 4 in cancer-endothelial co-culturing microenvironments promote proliferation, migration, and invasion of renal cancer cells. Cancer Cell Int. 2019;19:170.
[8]. Li M, Long X, Wan H, et al. Monocarboxylate transporter 1 promotes proliferation and invasion of renal cancer cells by mediating acetate transport. Cell Biol Int. 2021;45:1278–87.
[9]. Avitabile M, Succoio M, Testori A, et al. Neural crest-derived tumor neuroblastoma and melanoma share 1p13.2 as susceptibility locus that shows a long-range interaction with the
SLC16A1 gene. Carcinogenesis. 2020;41:284–95.
[10]. Song M, Zhong A, Yang J, et al. Large-scale analyses identify a cluster of novel long noncoding RNAs as potential competitive endogenous RNAs in progression of hepatocellular carcinoma. Aging (Albany NY). 2019;11:10422–53.
[11]. Lin HH, Tsai WC, Tsai CK, et al. Overexpression of cell-surface marker
SLC16A1 shortened survival in human high-grade gliomas. J Mol Neurosci. 2021;71:1614–21.
[12]. Jin Y, Qin X. Comprehensive analysis of transcriptome data for identifying biomarkers and therapeutic targets in head and neck squamous cell carcinoma. Ann Transl Med. 2020;8:282.
[13]. Al-Mustanjid M, Mahmud SMH, Royel MRI, et al. Detection of molecular signatures and pathways shared in inflammatory bowel disease and colorectal cancer: a bioinformatics and systems biology approach. Genomics. 2020;112:3416–26.
[14]. Eilertsen M, Andersen S, Al-Saad S, et al. Monocarboxylate transporters 1-4 in NSCLC:
SLC16A1 is an independent prognostic marker for survival. PLoS One. 2014;9:e105038.
[15]. Alexander J, Cukierman E. Cancer associated fibroblast: mediators of tumorigenesis. Matrix Biol. 2020;91-92:19–34.
[16]. Salimifard S, Masjedi A, Hojjat-Farsangi M, et al. Cancer associated fibroblasts as novel promising therapeutic targets in breast cancer. Pathol Res Pract. 2020;216:152915.
[17]. Payen VL, Mina E, Van Hee VF, et al. Monocarboxylate transporters in cancer. Mol Metab. 2020;33:48–66.
[18]. Philp NJ, Wang D, Yoon H, et al. Polarized expression of monocarboxylate transporters in human retinal pigment epithelium and ARPE-19 cells. Invest Ophthalmol Vis Sci. 2003;44:1716–21.
[19]. Demirbilek H, Hussain K. Congenital hyperinsulinism: diagnosis and treatment update. J Clin Res Pediatr Endocrinol. 2017;9(Suppl 2):69–87.
[20]. Kubelt C, Peters S, Ahmeti H, et al. Intratumoral distribution of lactate and the monocarboxylate transporters 1 and 4 in human glioblastoma multiforme and their relationships to tumor progression-associated markers. Int J Mol Sci. 2020;21:282.
[21]. Yu S, Wu Y, Li C, et al. Comprehensive analysis of the SLC16A gene family in pancreatic cancer via integrated bioinformatics. Sci Rep. 2020;10:7315.
[22]. Fei F, Guo X, Chen Y, et al. Polymorphisms of monocarboxylate transporter genes are associated with clinical outcomes in patients with colorectal cancer. J Cancer Res Clin Oncol. 2015;141:1095–102.
[23]. Xie J, Zhu Z, Cao Y, et al. Solute carrier transporter superfamily member
SLC16A1 is a potential prognostic biomarker and associated with immune infiltration in skin cutaneous melanoma. Channels (Austin). 2021;15:483–95.
[24]. Khan A, Valli E, Lam H, et al. Targeting metabolic activity in high-risk neuroblastoma through Monocarboxylate Transporter 1 (
SLC16A1) inhibition. Oncogene. 2020;39:3555–70.
[25]. Birsoy K, Wang T, Possemato R, et al.
SLC16A1-mediated transport of a toxic molecule is an effective strategy for targeting glycolytic tumors. Nat Genet. 2013;45:104–8.
[26]. Gray AL, Coleman DT, Shi R, et al. Monocarboxylate transporter 1 contributes to growth factor-induced tumor cell migration independent of transporter activity. Oncotarget. 2016;7:32695–706.
[27]. Chen X, Chen X, Liu F, et al. Monocarboxylate transporter 1 is an independent prognostic factor in esophageal squamous cell carcinoma. Oncol Rep. 2019;41:2529–39.
