Multi-Omics Analysis Identified TMED2 as a Shared Potential Biomarker in Six Subtypes of Human Cancer

Introduction

Among malignancies, cancer has the 2nd highest morbidity and mortality rates after cardiovascular diseases worldwide.¹ Metastasis is main characteristic of the malignant tumors which leads to poor prognosis of the cancer patients.² Currently, chemotherapy is the most widely used treatment approach for cancer patients; however, its efficacy is not up to the mark. On the other hand, immunotherapy is another kind of partially successful cancer treatment method, at present, it is very expensive,³ thus, targeted therapy remains the best treatment option for most cancer patients. Therefore, it is urgent for us to discover the underlying mechanisms of higher cancer incidence and detect some novel shared diagnostic and prognostic biomarkers that could aid in the targeted therapy of cancer patients.

Transmembrane p24 trafficking protein 2 (TMED2) is a member of the transmembrane emp24 domain (Tmed)/p24 family of proteins that are required for the trafficking of proteins among Golgi apparatus and the endoplasmic reticulum (ER).⁴ TMED2 normally functions during placenta and embryo morphogenesis,⁵ however, its role in the development and progression of tumorigenesis has been explored rarely. Few earlier studies have revealed the potential role of the TMED2 gene in different cancer subtypes including prostate cancer, colon cancer, and hepatocellular carcinoma.^6–8 However, best to our knowledge, the role of TMED2 in different other subtypes of human cancer is yet to uncover.

Recently, the development of advanced cancer expression databases and bioinformatics tools have helped cancer researchers to manipulate and analyze the experimentally obtained multi-omics data. To better understand the TMED2 role in the development and progression of human cancers, we analyzed its expression level variations, prognostic significance, promoter methylation level, genetic alterations, copy number variations (CNVs), pathways, CD8+ T immune cells infiltration, and gene-drug network using a comprehensive in silico approach. The findings of our study have provided some useful information regarding the correlation between TMED2 expression and its prognostic values in Cervical Squamous Cell Carcinoma (CESC), Esophageal Adenocarcinoma (ESCA), Head and Neck Cancer (HNSC), Kidney Renal Clear Cell Carcinoma (KIRC), Liver Hepatocellular Carcinoma (LIHC), and Lung Adenocarcinoma (LUAD), as well as suggested the role of TMED2 as a potential diagnostic and prognostic biomarker in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD.

Methods

Ethics Statement

This study was reviewed and approved by the ethical committee of the Islamia University of Bahawalpur, Bahawalpur, Pakistan.

UALCAN

Transcription expression level variations of TMED2 in different human cancers relative to normal controls were analyzed using UALCAN database (http://ualcan.path.uab.edu). This database is a reliable and easy-to-use platform for performing the multi-omics analysis of TCGA cancer data.⁹ In UALCAN, the transcription expression level of TMED2 was quantified as transcript per million (TPM) reads, and a Student’s t-test was employed for statistics purpose. A p-value (<0.05) was considered as statistically significant.

Kaplan–Meier Plotter

The Kaplan–Meier plotter (http://kmplot.com//analysis) tool has the capacity to compute prognostic values of more than 54K genes including mRNAs, miRNAs, and proteins in 21 subtypes of human cancer.⁹ In this study, the primary goal of utilizing this tool was to examine the clinical significance (overall survival [OS] and relapse-free survival [RFS]) of TMED2 in distinct subtypes of human cancer based on GEO, and TCGA datasets. A p-value, 95% confidence interval (CI), and hazard ratio (HR) were also determined and displayed.

GENT2 Database

GENT2 (http://gent2.appex.kr/) web server is developed for the multi-omics analysis of the cancer data acquired from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database.¹⁰ In our study, we utilized GENT2 for the validation of TMED2 expression using new independent cohorts of distinct cancer patients. For statistics, this tool uses a Student’s t-test, and a p-value of <0.05 was regarded as statistically significant.

Data Mining Through Human Protein Atlas Database

HPA (https:/www.proteinatlas.org/) is a database that aims to quantify and locate all the human proteins in cells, tissues, and organs using the immunohistochemistry (IHC) technique.¹¹ In the current study, HPA was utilized to examine the proteomics expression level of TMED2 in normal tissues and distinct subtypes of cancerous tissues. The observed protein expression level was graded as not detected, low, medium, and high, based on the intensity of staining and fraction of the stained cells.

MEXPRESS Database

MEXPRESS database was utilized to assess the Pearson correlation between TMED2 transcription expression and its promoter methylation levels in 309, 194, 542, 343, 414, and 477 samples of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD, respectively.¹² A p-value <0.05 was considered as significant.

The cBioportal Database

In this study, we utilized the cBioPortal database (http://cbioportal.org)¹³ to examine the TMED2-related genetic alterations and copy number variations (CNVs) across distinct human cancer subtypes using TCGA datasets. The OncoPrint option on the cBioPortal provides an overview of genetic mutations and CNVs in the gene(s) of interest above every sample bar.

TMED2 Protein–Protein Interaction Network Development and Pathway Analysis

A PPI network of TMED2-enriched genes was developed via STRING (https://string-db.org/)¹⁴ and visualized through Cytoscape software.¹⁵ Then, the pathway analysis of the TMED2-enriched genes was carried out using DAVID (http://david.ncifcrf.gov/summary.jsp)¹⁶ and a p value < 0.05 was considered as significant.

TMED2 Expression and Infiltrating Level of CD8+ T Immune Cells in Cancer Patients of Distinct Subtypes

TIMER (http://timer.cistrome.org/)¹⁷ is a database that aims to systematically analyze the immune infiltrates in distinct cancer subtypes (www.cistrome.shinyapps.io/timer) (48). TIMER applied the Spearman correlation method to infer the abundance of immune cells in association with gene expression using TCGA datasets. In distinct human cancer subtypes, the Spearman correlation between the TMED2 expression and CD8+ T immune cells was calculated in this study using TIMER. A p-value <0.05 was considered as significant.

Gene–Drug Interaction Network Analysis of TMED2

TMED2 gene–drug interaction network was constructed using the Comparative Toxicogenomics Database (CTD)¹⁸ and Cytoscape software. This network includes the information of chemotherapeutic drugs that could reduce or enhance the transcription or translation expression levels of TMED2.

Results

TMED2 Transcription Expression Analysis in Human Cancers and Normal Tissues

The transcription level of TMED2 in 24 major human cancer subtypes was explored via UALCAN database using TCGA datasets. Relative to solid normal tissues, TMED2 transcription level was significantly (p < 0.05) enhanced in all the 24 types of primary cancerous solid tissues (Figure 1) including Cervical Squamous Cell Carcinoma (CESC), Esophageal Cancer (ESCA), Head and Neck Squamous Cell Carcinoma (HNSC), Kidney Renal Clear Cell Carcinoma (KIRC), Liver Hepatocellular Carcinoma (LIHC), and Lung adenocarcinoma (LUAD) (Figure 1).

Figure 1 The analysis of difference in transcriptional level of TMED2 in different types of cancers tissue paired with normal controls. (A) Transcriptional level analysis of TMED2 across cancerous samples paired with normal controls, and (B) Transcriptional level analysis of TMED2 across cancerous samples only. Blue color represents the normal samples while red color indicates the cancer samples. *p < 0.05.

The Prognostic Potential of TMED2 in Various Types of Cancers

We further explored the correlation between OS, RFS and TMED2 expression in 24 major types of cancer patients using KM plotter. Elevated TMED2 expression was found to be significantly (p < 0.05) correlated with the reduced OS and RFS duration of the CESC (HR = 2.07, 95% CI: 1.29–3.32, P = 0.0021, HR = 1.79, 95% CI: 10.81–3.95, P = 0.014), ESCA (HR = 3.36, 95% CI: 1.68–6.73, P = 3e-04, HR = 1.6, 95% CI: 0.0–3.19, P = 0.018), HNSC (HR = 2.08, 95% CI: 1.53–2.83, P = 2.2e-06, HR = 2.75, 95% CI: 1.24–6.1, P = 0.0094), KIRC (HR = 1.64, 95% CI: 1.2–2.24, P = 0.0019, HR = 274586506.99, 95% CI: 0–Inf, P = 0.0021), LIHC (HR = 1.69, 95% CI: 1.18–2.43, P = 0.0038, HR = 1.36, 95% CI: 0.95–1.95, P = 0.05), and LUAD (HR = 1.65, 95% CI: 1.19–2.3, P = 0.0027, HR = 1.64, 95% CI: 0.95–2.83, P = 0.057) patients (Figure 2). Furthermore, the TMED2 dysregulation was also found associated with the OS and RFS of patients with distinct other types of cancers but the correlations were insignificant (p>0.05). Taken together, these data suggested that TMED2 might have a significant contribution to the development and progression of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD, thus the next part of our study will mainly focus on the unique role of TMED2 in these six types of human cancers.

Figure 2 Correlational analysis of TMED2 dysregulation with OS and RFS in distinct cancer subtypes. (A) Correlational analysis of TMED2 dysregulation with OS in distinct cancer subtypes, and (B) Correlational analysis of TMED2 dysregulation with RFS in distinct cancer subtypes. A p-value < 0.05 was considered as significant.

Association Between TMED2 Expression Level and Different Clinicopathological Features in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD

We further analyzed the TMED2 transcription expression levels in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients with respect to different clinicopathological features. We noticed that TMED2 transcription expression was also significantly up-regulated in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients relative to normal controls irrespective of different clinicopathological features including different cancer stages, races, genders, and age groups (Table 1). Clinicopathological distribution of the CESC, ESCA, HNSC, KIRC, LIHC, and LUAD cohorts are given in Tables 2–4.

Table 1 Clinicopathalogical Features-Specific Expression Pattern of TMED2 in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD Patients

Table 2 Clinicopathological Features of the CESC and ESCA Cohorts Included in the Present Study

Table 3 Clinicopathological Features of the HNSC and KIRC Cohorts Included in the Present Study

Table 4 Clinicopathological Features of the LIHC and LUAD Cohorts Included in the Present Study

TMED2 Transcription Level Validation Using New Independent Cohorts of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD

TMED2 transcription expression was re-analyzed using new expression data obtained through Affymetrix U133A and U133Plus2 microarray platforms from independent cohorts of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD via GENT2 database. Our results further validated the significant (p>0.05) up-regulation of TMED2 in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients relative to normal controls (Figure 3). Details of the new CESC, ESCA, HNSC, KIRC, LIHC, and LUAD datasets that were used for TMED2 transcription expression validation can be seen in Table 5.

Table 5 Details of the CESC, ESCA, HNSC, KIRC, LIHC, and LUAD Datasets Used for the TMED2 Expression Validation via GENT2 Web Server

Figure 3 Transcription expression level validation of TMED2 using independent CESC, ESCA, HNSC, KIRC, LIHC, and LUAD cohorts via GENT2 database. Blue color represents the normal samples while red color indicates the cancer samples. *p < 0.05.

Protein Expression Level Analysis of TMED2 in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD

TMED2 protein expression levels in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD and normal tissues were analyzed via the HPA database. As shown in Figure 4, an enhanced TMED2 protein expression was found in CESC (high), ESCA (high), HNSC (high), KIRC (high), LIHC (high), and LUAD (high) specimens relative to the paired normal CESC (medium), ESCA (medium), HNSC (medium), KIRC (low), LIHC (medium), and LUAD (medium). Taken together, our results have shown that TMED2 was also overexpressed at the translational level in cancer tissues of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD relative to the normal controls.

Figure 4 TMED2 translational expression analysis across distinct cancer tissues paired with normal controls via Human Protein Atlas (HPA) database (×200). (A) Cervical cancer, (B) Esophageal cancer, (C) Head and neck cancer, (D) Kidney cancer, (E) Liver cancer, and (F) Lung cancer.

Promoter Methylation Analysis of the TMED2 in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD

Promoter methylation performs the critical function in regulating gene expression, and any abnormality in promoter methylation can cause various genetic abnormalities including cancer.¹⁹ Herein, we utilized the MEXPRESS webserver to analyze whether TMED2 up-regulation is associated with its promoter methylation level or not. As per the analysis, we found a significant (p > 0.05) negative correlation between the promoter methylation and up-regulation of TMED2 in CESC, ESCA, HNSC, and KIRC, while a significant (p > 0.05) positive correlation in LIHC and LUAD. Collectively, these data suggested that TMED2 promoter hypomethylation might involve in its up-regulation in CESC, ESCA, HNSC, and KIRC (Figure 5) while hypermethylation of TMED2 in LIHC and LUAD (Figure 6) challenges the classical view of methylation and demand further extensive work for conclusive remarks.

Figure 5 A MRXPRESS based correlation analysis between TMED2 expression and its promoter methylation in CESC, ESCA, HNSC, and KIRC. (A) In CESC, (B) In ESCA, (C) In HNSC, and (D) In KIRC. A negative sign indicates the negative correlation between TMED2 expression and its promoter methylation using a specific probe at a specific CpG island while a positive sign indicates the positive correlation between TMED2 expression and its promoter methylation using a specific probe at a specific CpG island. *p < 0.05.

Figure 6 A MRXPRESS based correlation analysis between TMED2 expression and its promoter methylation in LIHC, and LUAD. (A) In LIHC, and (B) In LUAD. A negative sign indicates the negative correlation between TMED2 expression and its promoter methylation using a specific probe at a specific CpG island while a positive sign indicates the positive correlation between TMED2 expression and its promoter methylation using a specific probe at a specific CpG island. *p < 0.05.

Genetic Alteration Analysis of TMED2 in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD

Next, we determined the TMED2 genetic alteration frequencies in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD samples through cBioPortal webserver using TCGA PanCancer Atlas datasets. Missense alterations were more likely to observe in CESC, LIHC, while ESCA, HNSC, KIRC, and LUAD were more likely to harbor deep amplification and deep deletion, respectively (Figure 7). Taken together, it is speculated that low percentages of observed genetic alterations in CESC (0.4%), ESCA (0.5%), HNSC (0.6%), KIRC (0.3%), LIHC (0.3%), and LUAD (0.8%) samples might have a least participation in the expression regulation of TMED2.

Figure 7 TMED2-related genetic alterations and CNVs in TCGA CESC, ESCA, HNSC, KIRC, LIHC, and LUAD datasets. (A) In CESC dataset, (B) In ESCA dataset, (C) In HNSC dataset, (D) In KIRC dataset, (E) In LIHC dataset and, (F) In LUAD dataset.

PPI Network Construction, Visualization, and Pathway Enrichment Analysis of TMED2

We utilized STRING and Cytoscape resources to identify the TMED2 enriched genes. Functional interaction network analysis showed that TMED2 was closely associated with 16 different other genes (Figure 8A). We next performed the pathway analysis of TMED2 enriched genes using the DAVID tool. The obtained results showed that few TMED2 associated genes were significantly (p > 0.05) involved in 3 diverse pathways including Protein processing in endoplasmic reticulum, Legionellosis, and Endocytosis (Figure 8 and Table 6).

Table 6 Detail of Kyoto Encyclopedia of Genes and Genomes Pathway Analysis of the TMED2 Enriched Genes

Figure 8 A TMED2-enriched genes PPI network and pathway analysis. (A) A PPI network of TMED2 enriched genes, (B) Pathways analysis of the TMED2 enriched genes.

TMED2 Overexpression and Infiltrating Level of CD8+ T Immune Cells in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD Patients

To further explore whether TMED2 exerts any potential biological roles in CD8+ T immune infiltration or not, we carried out an integrated analysis using the TIMER database to find the association between TMED2 overexpression and CD8+ T immune cell infiltration in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD. As shown in Figure 9, a high level of TMED2 have a significant (p>0.05) negative correlation with CD8+ T immune cells level in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD (Figure 9).

Figure 9 Spearman correlational analysis of TMED2 expression with CD8+ T immune cells level in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD via TIMER. A p-value < 0.05 was considered as significant.

Gene–Drug Interaction Network Analysis of the TMED2

To find the association of TMED2 expression with different available cancer therapeutic drugs, we developed and analyzed the gene–drug interaction network via CTD database and Cytoscape software. In view of our results, it was observed that TMED2 expression can potentially regulate by a variety of drugs. For example, estradiol and bisphenol A can elevate the expression level of TMED2 while ketolides and glyphosate can reduce TMED2 expression level (Figure 10).

Figure 10 A gene–drug interaction network of TMED2 targeted drugs from CTD. Red arrows show chemotherapeutic drugs that can increase TMED2 expression while green arrows highlight chemotherapeutic drugs that can decrease TMED2 expression. The numbers of arrows between each interaction represent the supported numbers of studies.

Discussion

The TMED gene family was discovered around twenty years ago, its functions and underlying mechanisms are still unclear. Few earlier studies have reported that dysregulation of TMED genes can cause many diseases ranging from cancer to Alzheimer’s.²⁰ Some studies also revealed that TMED3 is a tumor suppressor gene and plays an important role in the initiation and development of colon cancer, prostate cancer, and hepatocellular carcinoma.^6–8 However, the role of TMED2 in cancer has been least studied so far. In the present study, we comprehensively analyzed TMED2 expression in 24 major subtypes of the human cancers to identify few TMED2-associated subtypes.

In our study, we observed that TMED2 was significantly (p > 0.05) overexpressed in all 24 major subtypes of human cancers relative to normal controls (Figure 1). Next, we also verified that TMED2 overexpression was significantly (p < 0.05) associated with the significant (p < 0.05) decreased in OS and RFS duration of the CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients. Taken together, these findings told us that TMED2 might play a crucial role in the initiation and progression of CESC, ESCA, HNSC, KIRC, LIHC, and LUAD, therefore, in the present study, our main focus are these six cancer subtypes. Next, we further re-examined TMED2 expression individually in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients of different clinicopathological features. A significant (p < 0.05) up-regulation of TMED2 was also observed in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients of various clinicopathological features in comparison to normal controls such as different cancer stages, patients races, patients genders, and patients ages.

The epigenome variations are thought to be linked with numerous types of diseases; however, DNA promoter methylation is the only epigenetic marker that is stable enough to detect in different kinds of diseased samples.²¹ Any abnormal change in the DNA promoter methylation level either hypo or hypermethylation might result in aberrant gene expression.²² Earlier studies have already identified the role of abnormal DNA promoter methylation in different types of human cancers.^23,24 In our study, we noticed that the TMED2 promoter methylation level was lower as per expectations in CESC, ESCA, HNSC, and KIRC samples relative to normal controls. In contrast, we observed that TMED2 promoter methylation level was unexpectedly higher in LIHC and LUAD samples as compared to the normal controls. Such differences in the promoter methylation levels of TMED2 may suggest that the TMED2 gene has different epigenetic regulations among different types of cancers. Taken together, we speculated that TMED2 promoter hypomethylation is closely linked with its overexpression in CESC, ESCA, HNSC, and KIRC, however, the causes of different TMED2 promoter hypermethylation in LIHC and LUAD are needed to be explored by further large-scale studies.

Genetic abnormalities are the molecular biomarker of tumorigenesis, and CNVs are the major types of genetic abnormalities that play a key role in the initiation and development of different human cancers.^25,26 In this study, we identified the minor percentages of the TMED2 genetic abnormalities and CNVs (0.4%, 0.5%, 0.6%, 0.%, 0.3%, and 0.8%) in the analyzed CESC, ESCA, HNSC, KIRC, LIHC, and LUAD samples, respectively. Overall, these results suggested that genetic abnormalities may not play any significant role in the expression regulation of TMED2, however, further detailed work is required to confirm it.

Although recent studies have identified several CESC, ESCA, HNSC, KIRC, LIHC, and LUAD associated diagnostic and prognostic biomarkers including different genes, such as RFC4, ADGRF4, ANXA8L1, TOP2A, HCAR3, MCM2, IRF6, and PDE2A in CESC,^27,28 COL1A2, DGCR8, COL1A1, POM121, TAF9, UPF3B, ZNF469, BCAP31, and COL3A1 in ESCA,^29,30 LAMC2, CCT3/4/5/6/7/8, and FGFR1-4 in HNSC,^31–33 PLCB2, TNFSF13B, VAV1, RAC2, and PARVG in KIRC,^34,35 CDC20, MTHFD1L, CDCA8, CDCA5, KIF2C, and KIFC1 in LIHC,^36,37 CCNB1, SPP1, TOP2A, MKI67, CHEK1, RRM2, and CDK1 in LUAD.³⁸ However, none of these or any other biomarkers have been generalized so far in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients of different clinicopathological features. In the current study, we have noticed the significant (p < 0.05) TMED2 up-regulation in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients of different clinicopathological features including different cancer stages, patients races, genders, and age groups relative to controls. Additionally, the promoter methylation level, OS, and RFS information of TMED2 have also highlighted its potential as a useful novel potential biomarker in these cancers. Best to our knowledge, this study is the first to report a shared clinicopathological features-specific diagnostic and prognostic role of TMED2 in six different cancers including CESC, ESCA, HNSC, KIRC, LIHC, and LUAD, which may help to design the effective diagnostic and therapeutic strategies for cancer patients.

Recent studies have investigated that some genes regulate the tumor microenvironment by changing the levels of immune cells.^39,40 Therefore, we analyzed the association between TMED2 and CD8+ T immune cells infiltration in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD. As a result, TMED2 overexpression has shown a negative correlation with CD8+ T immune cells level in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD. Collectively, this correlation scenario between TMED2 overexpression and CD8+ T immune infiltration indicates that TMED2 may initiate tumor development by changing the level of immune infiltration in the tumor microenvironment, but further work is required to confirm it.

The PPI network of TMED revealed that it directly interacts with the 16 different genes and a few of them were involved in 3 important signaling pathways including Protein processing in endoplasmic reticulum, Legionellosis, and Endocytosis (Table 6). Additionally, we have also identified few potential drugs that could regulate TMED2 expression, suggesting their importance as effective participants in designing the most appropriate therapeutic strategies for CESC, ESCA, HNSC, KIRC, LIHC, and LUAD (Figure 10).

Conclusion

In conclusion, we found that TMED2 is overexpressed in 24 major types of human cancers while it significantly reduced the OS and RFS duration of only CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients. This implies TMED2 is involved in the oncogenesis of these cancers. In addition, further multi-omics analysis also suggested TMED2 as the novel diagnostic and prognostic biomarker of survival in CESC, ESCA, HNSC, KIRC, LIHC, and LUAD patients of different clinicopathological features. However, additional work needed to be done on a large scale prior to clinical implications.