1. TMEM33 is highly expressed in cervical squamous cell carcinoma and endocervical adenocarcinoma
Based on the TCGA and GTEx databases, we examined the expression profile of TMEM33 in pan-cancer analysis. As shown in Figure 1a, among 33 cancer types, TMEM33 was significantly highly expressed in 24 cancers compared with normal tissues. For TCGA tumors and cancer-adjacent normal tissues, TMEM33 expression was significantly up-regulated in 14 cancer types (Fig. 1b). Specifically, TMEM33 expression was increased in CESC (Fig. 1c), and was further validated in the GSE63514 dataset (Fig. 1d). The expression of TMEM33 was higher in adenosquamous compared with squamous cell carcinoma of CESC (Fig. 1e). In order to explore the clinical value of TMEM33 in CESC diagnosis, we performed ROC curve analysis. As the area under the curve (AUC) was 0.881, TMEM33 showed significantly high sensitivity and specificity for CESC diagnosis (Fig. 1f).
2. Association of TMEM33 expression and clinicopathological characteristics in CESC
We examined the clinicopathological characteristics of CESC patients based on the mean value of TMEM33 expression. As shown in Table 1, the high expression of TMEM33 was significantly associated with Height (P=0.048), N stage (P=0.003) and Histological type (P=0.003), while it was not associated with other features. The results of univariate analysis using logistic regression indicated that TMEM33 expression was correlated with poor prognostic clinical characteristics in CESC patients (Table 2). High TMEM33 expression was correlated with N stage (OR=0.372, P=0.002) and Histological type (OR=0.366, P=0.002).
3. Identification of TMEM33 as an independent prognostic indicator in CESC
Next,Kaplan-Meier analysis was used to assess the predictive value of TMEM33 on clinical outcomes. As shown in Figures 2a-c, high expression of TMEM33 predicted poor prognosis in both overall survival (OS, HR: 1.97, P=0.01), progress free interval (PFI, HR: 3.13, P=0.001) and disease specific survival (DSS, HR: 2.40, P=0.006). In addition, we conducted univariate and multivariate analysis using the Cox proportional hazards regression model. As shown in Table 3 and Supplementary Figure 1, the univariate cox regression analysis revealed TMEM33, clinical T, N and M stages were positive prognostic factors informing worsen CESC patient overall survival. However, only T stage, N stage and TMEM33 were independent risk factors for overall survival in multivariate regression cox analysis.
Moreover, a nomogram was constructed taking into account the statistically significant prognostic factors in multivariate Cox regression analysis to predict 3, 5 and 7-year overall survival (Fig. 2d). The calibration plots demonstrated good agreement between prediction and observed events (Fig. 2e). These data illustrate that the high expression of TMEM33 predicts poor prognosis of CESC.
4. Functional enrichment analysis of differentially expressed genes of TMEM33
To better understand the potential function and molecular pathways of the TMEM33 gene in CESC, we analyzed the differentially expressed genes (DEGs) based on the TCGA dataset. A total of 1,178 genes related to TMEM33 expression were altered (711 up-regulated, 467 down-regulated, |log2(FC)|>1 & p.adj<0.05), which may reflect the potential value of TMEM33 on CESC pathogenesis (data not shown).
Subsequently, TMEM33 and its DEGs were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. The top 10 enrichment pathways from each analysis result were exhibited. As shown in Figure 3a, the gene clusters were involved in epidermal cell differentiation, production of the molecular mediator of immune response, immunoglobulin production, and complement activation in the category of biological processes. For molecular function, the DEGs were enriched at receptor ligand activity, glycosaminoglycan binding, immunoglobulin receptor binding and antigen binding, etc. (Fig. 3b). In the category of cellular components, the DEGs were mostly distributed in immunoglobulin complex, apical part of cell, apical plasma membrane, intermediate filament cytoskeleton and blood microparticle (Fig. 3c). The KEGG pathway analysis showed that the DEGs were enriched in neuroactive ligand-receptor interaction, estrogen signaling pathway, vascular smooth muscle contraction, etc. (Fig. 3d). The GO and KEGG analysis demonstrate that TMEM33 might participate in signal transduction pathways and immune regulation.
5. GSEA identifies a TMEM33-related signaling pathway
Using the low- and high-TMEM33 expression datasets, we performed GESA to identify the most significantly altered pathways based on the gene sets from the MSigDB (c2.cp.kegg.v6.2.symbols.gmt). As shown in Figure 4a-f, gene sets related to genes encoding extracellular matrix proteins, cellular response to external stimulus, calcium mobilization, innate immune system, adaptive immune system and complement cascade, with the threshold of NES >1.5 and p.adj<0.05. Interestingly, in addition to the pathways that have been reported to reflect the biological functions of TMEM33 as a transmembrane protein, some immune-related signaling pathways were enriched, which provided additional evidence that TMEM33 might be implicated in immune response.
6. Association of TMEM33 and immune cell infiltration in CESC
To confirm our hypothesis, we determined the infiltration of 24 immune cell types in CESC using the ssGSEA method, and subsequently the association between TMEM33 and immune cell infiltration was investigated by Spearman’s correlation analysis. As shown in Figure 5a, T helper cells and Eosinophils were positively correlated with TMEM33 expression. However, Th1 cells, Treg, iDC, DC, Cytotoxic cells, B cells, T cells, CD56dim NK cells, pDC and aDC showed a significant negative association with TMEM33. More specifically, we assessed the infiltration levels of the top two relevant immune cells of each group, T helper cells, Eosinophils, Th1 cells and Treg in distinct TICRR groups, which showed consistent trends with the forest plots (Fig. 5b). Further, the correlation between TMEM33 expression and T helper cells, Th1 cells, Treg, DC, Cytotoxic cells and B cells in CESC were exhibited by scatter plots (Fig. 5c).
7. Functional inference of TMEM33 in CESC
In order to investigate the internal mechanism of TMEM33 in tumorigenesis, the PPI network analysis was performed by utilizing the STRING database. As shown in Figure 6a, the interaction network of 50 TMEM33-binding proteins with the experimental evidence identification was visualized. In addition, we analyzed the similar genes of TMEM33 in CESC using the “similar gene detection” module of GEPIA2. Corresponding hierarchical clustering analysis of the top 20 similar genes was displayed by the heatmap (Fig. 6b). Through comparing the 50 TMEM33-interacted genes and the 200 TMEM33-similar genes, we screened out the common molecule RNF4 (Fig. 6c), which plays a role in protein ubiquitination[17]. Moreover, the TMEM33 expression level was remarkably positively correlated with similar genes of RNF4, OCIAD1, TMED5, DHX15, MED28 and LETM1 (Fig. 6d), which have been reported to be implicated in tumorigenesis of various cancers[18-23].
We further performed enrichment analysis using the above 250 TMEM33-interacted and similar genes. The top 5 enrichment pathways of MF, CC and BP were shown in Figure 6e. However, there was no pathway enriched with KEGG owing to the few number of genes selected.
8. TMEM33 promotes cell proliferation of cervical cancer cells in vitro
To investigate the expression level of TMEM33 in cervical cancer cells, we performed RT-qPCR and immunoblotting in cervical cancer cell lines HeLa, SiHa, CaSki, H8 and C33A. TMEM33 was up-regulated in SiHa at the mRNA level (Fig. 7a), while highly expressed in HeLa and SiHa at the protein level (Fig. 7b). Based on that analysis, HeLa (HPV 18+) and SiHa (HPV 16+) cells were chosen for further functional studies. TMEM33 knockdown was performed via RNAi and the transfection efficiency was shown in Figure 7c.
Since six tumorigenesis related TMEM33-correlated genes were filtered by bioinformatics tools (Fig. 6c), we verified their association with TMEM33 in cervical cancer cells. As shown in Figure 7d and e, the mRNA expression of RNF4, OCIAD1, TMED5, MED28 and LETM1 were significantly decreased in HeLa-siTMEM33 cells compared with control cells, while all six genes were significantly down-regulated in SiHa-siTMEM33 cells, indicating a certain correlation existed.
Next, cell viability was assessed by CCK-8 assay. HeLa-siTMEM33 and SiHa-siTMEM33 cells showed a significant reduction of cell viability compared with control cells after transfection in 96 hours (Fig. 7f-g). In addition, colony formation assay showed that TMEM33 knockdown reduced the clonogenicity of HeLa cells (Fig. 7h) and SiHa cells (Fig. 7i). The above results suggest that TMEM33 contributes to cell proliferation in cervical cancer HeLa and SiHa cells, which provide evidence of the carcinogenesis role of TMEM33 in CESC.