Multiomics characteristics and immunotherapeutic potential of EZH2 in pan-cancer

With more than 19 million new cases and nearly 1000000 deaths, malignancies continue to pose a severe threat to public health worldwide [1]. Despite substantial advances in surgical techniques, drug treatments, and screening means for several cancers, the dismal clinical outcome of massive cancer patients remains unchanged [2–5]. Given the complexity of cancers and the holistic nature of the human body, a comprehensive pan-cancer analysis of some critical oncogenes is of great importance to investigate the molecular mechanisms of cancer development and promotion, thereby improving the prognosis of patients. The Cancer Genome Atlas (TCGA) and Gene Expression Ontology (GEO) projects include bulk genomic and clinical datasets of various tumors, providing a solid foundation for our subsequent pan-cancer analysis [6,7].

Enhancer of zeste homolog 2 (EZH2) belongs to the polycomb group genes (PcGs), and it is a class of essential epigenetic regulator [8]. Generally, EZH2 can regulate the expression of downstream target genes via the following three mechanisms, including polycomb repressive complex 2 (PRC2)-dependent trimethylation of Lys-27 in histone 3 (H3K27me3), PRC2-independent methylation, and PRC2- and methylation-independent gene transactivation [9–15]. When EZH2 is mutated or aberrantly expressed, it usually indicates the development and progression of cancers. The abnormal overexpression of EZH2 was first observed in prostate cancer and correlated with shorter survival time [16]. Afterward, up-regulation of EZH2 was also detected in endometrial carcinoma, anaplastic thyroid carcinoma, esophageal cancer, nasopharyngeal carcinoma, breast cancer, and gastric cancer (GC) [17–23]. As a key tumor suppressor of P21, EZH2 can silence the transcription of P21 by combining with its promoter region, leading to the proliferation of GC cells [24]. EZH2 could induce the oncogenic transformation of mammary epithelial cells and serve as a key biomarker of aggressive breast cancer [25]. By enhancing the expression of CCL5, EZH2 facilitated the recruitment and aggregation of macrophages in lung cancer, led to the invasion of tumor cells, and accelerated the epithelial–mesenchymal transition (EMT) process of pancreatic cancer cells [26,27]. In addition, increased expression of endothelial EZH2 could promote tumor angiogenesis through methylation and silencing of vasohibin1 [28]. Given the essential role of EZH2 in the development and progression of cancer, therapeutic approaches targeting EZH2 have emerged as a promising approach for treating multiple cancer, such as inhibiting the activity of EZH2 methyltransferase, degradation of EZH2, breaking the structure of PRC2, and combining EZH2 inhibitors with other anticancer regimens [29–36]. All the studies above reveal that EZH2 may serve as a potential pan-cancer marker for prognosis prediction and improved treatment.

In the present study, the structure, expression, genetic mutation, DNA methylation, protein phosphorylation, and prognostic significance of EZH2 in pan-cancer were investigated to elucidate the potential oncogenic mechanism in various cancers. Moreover, we also explored and evaluated the correlation of EZH2 with immune escape and several important immunotherapy response biomarkers.

The chromosome localization information of EZH2 was retrieved and viewed through the GeneCards website (https://www.genecards.org/). The EZH2 gene sequence (including coding DNA sequences, UTR regions, exons, and introns) was downloaded via the National Center for Biotechnology Information (NCBI) Gene database (https://www.ncbi.nlm.nih.gov/gene), and its protein composition structure (involving domains, region, compositional bias) was searched from the UniProt database (https://www.uniprot.org/). Afterward, the illustrator of biological sequences (IBS) online tool (http://ibs.biocuckoo.org/) was applied to visualize amino acid sequences and nucleotides of EZH2 [37]. Additionally, the NCBI HomoloGene database (https://www.ncbi.nlm.nih.gov/homologene) was used to investigate the evolution of the EZH2 family, and a phylogenetic tree based on EZH2 protein was created to elucidate the evolutionary associations among different biological species. The intracellular distribution of EZH2 protein was obtained by retrieving the ‘SUBCELL’ module of the Human Protein Atlas database (HPA, https://www.proteinatlas.org/).

The mRNA expression levels of EZH2 in normal human tissues were investigated via the Genotype-Tissue Expression (GTEx) portal (https://commonfund.nih.gov/GTEx/), and its expression distribution in malignant tumor tissues was detected using the HPA database. Four types of organs with the most significant mRNA expression of EZH2 in normal and tumor tissues were screened out. Immunohistochemistry (IHC) images from the HPA website were extracted to further demonstrate the expression levels of EZH2 in these organs. Moreover, the expression distribution of EZH2 in normal human tissues and cancerous cell lines was also obtained and explored through the HPA web server. Based on the GEPIA2 online tool (http://http://gepia2.cancer-pku.cn), the expression of EZH2 in tumors and matched paracancerous tissues at the mRNA level was compared [38]. UALCAN is a comprehensive, publicly convenient, and interactive portal that facilitates deep analyses of cancer omics data, including the clinical proteomic tumor analysis consortium (CPTAC) data (http://ualcan.path.uab.edu/) [39]. The differential expression of EZH2 in primary cancer and corresponding adjacent normal tissues at the protein levels was examined based on the CPTAC dataset.

Information on EZH2 genetic alteration features of all TCGA pan-cancer samples was queried by logging into the cBioPortal database (http://www.cbioportal.org/), such as mutation frequencies, alteration types, copy number alteration (CNA), and the mutated site information of EZH2 was marked on the whole amino acid sequence and three-dimensional (3D) structure. The catalogue of somatic mutations in cancer (COSMIC) is an online website resource involving the bulk somatic mutation data in various cancers, which was exploited to acquire detailed single nucleotide variants (SNV) materials of EZH2 (http://cancer.sanger.ac.uk/cancergenome/projects/cosmic/) [40]. The mutational landscape of EZH2 in all TCGA pan-cancer samples was summarized and delineated via Sangerbox 3.0 platform (http://sangerbox.com/). To assess the effects of CNA and SNV events on EZH2 expression, the association between mRNA expression levels of EZH2 and CNV/SNV was explored through Pearson correlation analysis based on the GSCA database (http://bioinfo.life.hust.edu.cn/GSCA/).

The methylation levels of the EZH2 promoter were compared between TCGA pan-cancer and relevant normal tissues by querying the ‘methylation’ module of the UALCAN web portal [41]. The promoter methylation level is expressed by β values, with 0.25–0.3 being considered hypermethylation and 0.5–0.7 being considered hypomethylation [41,42]. The online platform MethSurv (https://biit.cs.ut.ee/methsurv/) was applied to explore the methylation levels of single CpG sites with EZH2 in 25 types of malignancies [43]. In addition, the Sangerbox 3.0 online server was used to assess the correlation between EZH2 expression and RNA methylation regulators such as N1-methyladenosine (m1A), 5-methylcytosine (m5C), and N6-methyladenosine (m6A).

The protein phosphorylation sites of EZH2 from amino acid sequences were searched and mapped by retrieving the PhosphoSitePlus network (http://www.phosphosite.org). The protein phosphorylation data of CPTAC samples in the UALCAN database were adopted to compare phosphorylation levels of EZH2 protein between tumor and matched paracancerous tissues.

The ‘Survival Map’ module from the GEPIA2 online tool was used to select all TCGA malignancies in which EZH2 expression significantly correlated with survival outcomes, including disease-free survival (DFS) and overall survival (OS). For these cancer types, samples were categorized into a high-expression or a low-expression subgroup for subsequent Kaplan–Meier (K–M) survival analysis according to EZH2 expression. Furthermore, a univariate regression analysis was performed to evaluate the effect of EZH2 expression on disease-specific survival (DSS) and progression-free survival (PFS) of patients in the TCGA pan-cancer cohort.

The expression levels of EZH2 in various pathological stages of TCGA pan-cancer samples were obtained and compared by querying the ‘Stage Plot’ module of GEPIA2. The link between the transcript expression levels of EZH2 and TMN staging was investigated and probed through the Sangerbox 3.0 website. To elucidate the relationship between EZH2 aberrant expression and tumor metastasis, the TNM plot web server (https://tnmplot.com/analysis/) was used to assess the expression of EZH2 in tumor-adjacent normal tissues, primary tumors, and metastases.

The TIMER algorithm in the Sangerbox 3.0 database was applied to examine the correlation between EZH2 expression and the infiltration abundance of the six most common immune cells, including CD4 + T cells, CD8 + T cells, B cells, macrophages, dendritic cells (DC), and neutrophils [44]. Based on Spearman correlation analysis, scatter plots were generated to depict the relationship between the expression of EZH2 and infiltrations of six immune cells. In addition, the immune estimation resource TIMER2 (http://timer.cistrome.org/) was adopted to assess the correlation of EZH2 expression with the infiltration levels of three immunosuppressive cells, namely cancer-associated fibroblasts (CAFs), regulatory T cells (Tregs), and myeloid-derived suppressor cells (MDSCs). To further demonstrate the link between EZH2 expression and the tumor-immunosuppressive microenvironment, the tumor immune dysfunction and exclusion (TIDE) website (http://tide.dfci.harvard.edu) was exploited to evaluate the effect of altered EZH2 expression on T-cell dysfunction phenotype [45].

Recently, immunotherapy has emerged as a powerful next-generation approach for treating various malignant tumors [46,47]. Several significant biomarkers, such as immune checkpoint molecules, tumor mutation burden (TMB), microsatellite instability (MSI), mismatch repair (MMR) status, and neoantigens, were closely related to immunotherapy-associated response [48–52]. To further explore the value of EZH2 in predicting immunotherapy response, the Spearman correlation coefficient was employed to assess the association between the EZH2 expression and these known predictive biomarkers for immunotherapy.

Based on the Genomics of Drug Sensitivity in Cancer (GDSC) pharmacogenomics database (https://www.cancerrxgene.org/), the sensitivity and resistance of anticancer agents targeting EZH2 were determined and examined. The half-maximal inhibitory concentration (IC50) values of the two most effective anticancer drugs targeting EZH2 were calculated and compared between EZH2 mutant and wild types. Moreover, the chemical formulas and molecular constructs of the two most efficient anticancer drugs were obtained from the DrugBank website (https://www.drugbank.ca/).

A protein–protein interaction (PPI) network was built through the online server GPS-Prot (http://gpsprot.org/index.php) to screen out the top 19 hub genes that interacted closely with EZH2 (confidence > 0.85) [53]. Subsequently, a total of 20 genes (including EZH2) were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and Gene Ontology (GO) functional annotation analyses. The GeneMANIA tool (http://genemania.org/) was employed to decipher the biological functions of EZH2 [54]. Based on the median expression levels of EZH2, all patients in the TCGA pan-cancer cohort were dichotomized into EZH2 low‐ or high‐expression subgroups for further gene set enrichment analysis (GSEA) [55].

The Starbase (http://starbase.sysu.edu.cn/index.php) and miRNet (https://www.mirnet.ca/) were used to evaluate and predict miRNAs targeting EZH2. The target miRNAs for EZH2 were identified by capturing the intersection of the results predicted by the above databases. With the help of TargetScanHuman (http://www.targetscan.org), the complementary sequences of EZH2 and target miRNAs were then revealed. Moreover, long noncoding RNAs (lncRNAs) targeting EZH2 miRNAs were determined by the miRNet database to create a competing endogenous RNA (ceRNA) regulatory network.

For the validation of mRNA levels, quantitative real-time polymerase chain reaction (qRT-PCR) assays were carried out in accordance with the manufacturer’s protocols. Normal digestive system epithelial cell and corresponding cancer cells were purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control. Forward and reverse primers applied for amplification were, respectively, EZH2-F (TGGTGCTGAAGCCTCAATGT), EZH2-R (CGGGAGCTGGAGCTATGATG); GAPDH-F (CCCACTCCTCCACCTTTGAC), GAPDH-R (CCACCACCCTGTTGCTGTAG). PCR data were processed using Microsoft Excel, and graphs were drawn using GraphPad Prism 7.0. The relative expression of mRNA was computed by using the 2^−ΔΔCq method.

To further confirm the differential expression of EZH2 at the translation levels, we collected four pairs of paraffin-embedded aerodigestive tract malignancies and matched paracancerous normal tissues from the First Affiliated Hospital of Nanchang University, namely LIHC, STAD, COAD, and PAAD. The paraffin-embedded tissues were cut into 4 µm sections, then deparaffinized for 2 h at 80°C and hydrated. These samples were incubated with 1:50 dilution of anti-EZH2 rabbit mAb (CST, #5246). After 1 h incubation with antirabbit secondary antibody at room temperature, the color of the antibody staining was revealed using diaminobenzidine (DAB). Finally, the stained slides were observed under a microscope at 200× magnification for pathologic assessment.

The present study aims to elucidate how EZH2 affects the development and progression of various human cancers. The flowchart presented the entire analysis process throughout the study (Figure 1A). The human EZH2 (Gene ID: 2146) gene is mapped on the seventh chromosome at region q36.1, including 20 exons and 19 introns with a length of 76978 bp (Figure 2A). This gene encodes a 746 amino acids protein, consisting of cysteine-rich (CXC) and SET domains, five regions, and four compositional biases (Figure 2B). The 3D structures of EZH2 protein were predicted by the homology database ModBase, as shown in Figure 2C. Next, the protein sequences coded by EZH2 mRNA from ten various species were compared to detect similar structures, such as WD repeat-binging protein EZH2 (ptfam11616) and SET domain (cl02566) (Figure 2D). Furthermore, a phylogenetic tree was established to further investigate the evolution characteristics of EZH2 protein and its homologs among different species (Figure 2E).

Additionally, EZH2 could encode five isoforms involving histone-lysine N-methyltransferase EZH2 isoforms a–e, predominantly located in the nucleoplasm. For example, within the nucleoplasm of REH (human acute lymphocyte leukemia cell), U-2 OS (human osteosarcoma cell), and SiHa (human cervical squamous cancer cell) (Figure 2F).

According to the GEPIA2 online tool, the ubiquitously spread of EZH2 in multiple human cancer and normal tissues can be observed (Figure 3A). To determine the specific mRNA expression of EZH2 in normal tissues, the mRNA expression data of EZH2 from the GTEx database were employed to examine the expression status of EZH2. The results revealed that the highest mRNA expression levels of EZH2 in normal human tissues were detected in the testis, followed by (in descending order) the esophagus, the skin, and the spleen (Figure 3B). In addition, IHC-staining images from the HPA database were conducted to further display the protein expression of EZH2 in the above four normal organs (Figure 3C). The expression levels of EZH2 in tumor tissues were also assessed through the HPA website. It can be seen that the EZH2 was most expressed in lung cancer, followed by colorectal cancer, head cancer, neck cancer, testis cancer, and renal cancer (Figure 4A). IHC-staining images were used to confirm strong positive protein expression of EZH2 in the top-ranked four human cancerous organs, with EZH2 expression being the most prominent (Figure 4B). Besides, the expression distribution of EZH2 in normal human tissues and cancerous cell lines was also evaluated. The results show that gastric mucous cells reached the highest EZH2 expression levels in normal human tissue cell lines (Figure 3D), while lymphoid cancer cells (e.g., MOLT-4, HDLM-2) were highest in cancer cell lines (Figure 4C).

Due to few normal tissues in the TCGA database, the mRNA expression data of normal human tissues derived from the GTEx database were integrated to test the differences in mRNA expression of EZH2 between tumor and matched paracancerous tissues. In addition to the dramatic reduction in acute myeloid leukemia (LAML), the mRNA expression levels of EZH2 were markedly elevated in various cancer specimens, including bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), glioblastoma multiforme (GBM), brain lower-grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung squamous cell carcinoma (LUSC), ovarian serous cystadenocarcinoma (OV), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), and uterine carcinosarcoma (UCS), compared with corresponding normal counterparts (Figure 5A). Compared with normal breast tissues, the protein expression levels of EZH2 were dramatically up-regulated in pancreatic adenocarcinoma (PAAD), head and neck squamous cell carcinoma (HNSC), lung adenocarcinoma (LUAD), UCEC, LIHC, and GBM, while the expression of EZH2 was significantly lower in breast cancer specimens (Figure 5B). Moreover, higher expression of EZH2 was detected in patients with MSI-High (MSI-H) than those with MSI-Low (MSI-L). In summary, EZH2 was abnormally overexpressed in multiple human cancer tissues, indicating that it may function as an oncogene in various tumor development and progression.

Genetic alterations are strongly related to the development of malignancy, and mutated genes can be potential targets for the therapy of cancers [56,57]. Given that EZH2 gene mutations might serve as a promising molecular target for the treatment of various human cancers, the mutational profile of EZH2 in human cancers was investigated using TCGA pan-cancer samples. Among 10953 TCGA pan-cancer patients, 290 (2.6%) experienced EZH2 mutation, with amplification and missense mutation being the predominant types of EZH2 gene alteration (Figure 6A). The highest mutation frequency of EZH2 (>8%) was discovered in TCGA UCEC samples, in which ‘mutation’ was the primary type, followed by OV (>7%) (Figure 6B). Notably, all uveal melanoma (UVM) cases with EZH2 gene alterations had only the ‘deep deletion’ type. The types and sites of EZH2 genetic mutations were noted throughout the whole amino acid sequence, and the E745k site with the highest alteration frequency was detected in four cases of UCEC, one case of STAD, and one case of colorectal adenocarcinoma (Figure 6D). As displayed in Figure 6C, the mutated site information of EZH2 was further exhibited by 3D structure. In addition, a more detailed mutational landscape of EZH2 in TCGA pan-cancer samples was explored and depicted through the Sangerbox 3.0 database (Supplementary Figure S1D).

Given that ‘missense mutation’ of SNV and ‘amplification’ of copy number variation (CNV) were the major types of EZH2 mutation, it was necessary to conduct further exploration and investigation on the SNV and CNV features of EZH2. Subsequent results suggested that missense mutation was the most common type of EZH2 SNV, and the highest proportion of SNV categories was G > A (19.90%), followed by C > T (15.52%), A > T (11.59%), and T > A (10.71%) (Supplementary Figure S1A,B). The SNV percentage of UCEC was the highest, accounting for 44%, followed by skin cutaneous melanoma (SKCM) (Supplementary Figure S1F). For CNV, diploid, gain, amplification, and shallow deletion were more prevalent, while deep deletion was rare (Supplementary Figure S1C). There was a strong positive association between EZH2 expression and CNV in most tumors (Supplementary Figure S1E). The CNV pie chart also revealed that heterozygous amplification was concentrated in various cancers, except for LAML (Supplementary Figure S1E). The above analysis indicated that different EZH2 mutational signatures in various cancers might be closely associated with the abnormal expression of EZH2.

DNA methylation patterns have been reported to contribute to tumorigenesis [58–60]. To examine whether DNA methylation of EZH2 was associated with tumorigenesis, the ‘TCGA methylation’ module of the UALCAN database was used to assess the promoter methylation level of EZH2. The results suggested that compared with corresponding normal tissues, the methylation state of EZH2 was lower in BLCA, UCEC, LUAD, LUSC, prostate adenocarcinoma (PRAD), READ, testicular germ cell tumors (TGCT), and thyroid carcinoma (THCA), and higher in cholangiocarcinoma (CHOL) (Figure 7A). Further analysis demonstrated that EZH2 expression was positively correlated with RNA methylation regulators (m1A, m5C, m6A) (Supplementary Figure S2). As shown in Supplementary Figure S3, the methylation levels of single CpG sites were correlated with EZH2 in 25 malignancies, and the highest DNA methylation level was easily observed at the cg18416251 site.

Protein phosphorylation is one of the most critical post-translational modifications that can significantly affect progression through multiple pathways, including MAPK, tyrosine kinases, PI3K/Akt, and TGF-β signaling [61]. To understand the significance of EZH2 protein phosphorylation in pan-cancer, the phosphorylation sites in EZH2 protein were mapped via the PhosphoSitePlus website. As depicted in Figure 7B, the most predominant phosphorylation locus for EZH2 protein was located at T487 (flanking sequence: APAEDVDtPPRKKKR). Afterward, the phosphorylation levels of EZH2 at the T487 site were examined and compared between the tumor and adjacent normal tissues. GBM (P=1.13 × 10^-2), LIHC (P=1.6 × 10^-18), LUAD (P=4.73 × 10^-24), and breast cancer (P=4.73 × 10^-8) had higher phosphorylation degree in T487 site of EZH2 protein, revealing that protein phosphorylation of EZH2 at T487 site might serve as a facilitator in the development and progression of these cancers (Figure 7C).

To explore the prognostic value of EZH2 in pan-cancer patients, the association between EZH2 expression levels and clinical outcome (OS and DFS) of cancer patients was evaluated by retrieving the ‘Survival Map’ module of the GEPIA2 database. As illustrated in Figure 8A, overexpressed EZH2 was associated with unfavorable OS of patients with adrenocortical carcinoma (ACC), kidney renal clear cell carcinoma (KIRC), LGG, LIHC, mesothelioma (MESO), pheochromocytoma and paraganglioma (PCPG), and PRAD. Besides, overexpressed EZH2 was a protective factor for THYM patients. High expression of EZH2 was also a significant risk factor for DFS in patients with ACC, BLCA, KICH, kidney renal papillary cell carcinoma (KIRP), LGG, LIHC, PRAD, and THCA (Figure 8B). Based on the intervention of noncancer-associated deaths, the correlation between the expression level of EZH2 and DSS in pan-cancer was examined using Cox regression analysis. The results showed that EZH2 expression affected DFI in patients with glioma (GBMLGG), pan-kidney cohort (KIPAN), LGG, ACC, KICH, LIHC, MESO, PCPG, KIRP, KIRC, PRAD, UVM, PAAD, and SKCM-P (Supplementary Figure S4B). Furthermore, whether the expression status of PCSK9 was associated with PFS in 38 cancer types was also investigated. Aberrant EZH2 expression was associated with PFS in patients with GBMLGG, ACC, KIPAN, PRAD, LIHC, LGG, KICH, UVM, THCA, KIRP, PCPG, MESO, and PAAD (Supplementary Figure S4A). The above results indicated that abnormal expression of EZH2 was closely related to the survival outcome of various cancers, especially ACC, LGG, PRAD, and LIHC.

To further exploit the potential value of EZH2 in guiding clinical decisions, the association between the expression levels of EZH2 and clinicopathological features was explored and assessed. By examining the expression levels of EZH2 in different pathological stages, it was observed that EZH2 expression was up-regulated in higher pathological stages of ACC, kidney chromophobe (KICH), and KIRC, while it was decreased in LIHC (Supplementary Figure S5A). EZH2 expression varied with different T stages of LAML, PAAD, and ACC, and the expression levels differed at various M stages of GBMLGG, HNSC, LIHC, and PAAD (Supplementary Figure S5C,D). Moreover, the expression levels of EZH2 were significantly increased in lung, liver, prostate, and kidney metastatic tissues compared with corresponding primary tumor tissues (Supplementary Figure S5B). These results demonstrated that EZH2 might function as one oncogene in the progression and metastasis of various tumors.

Innate immune cells (e.g., macrophages, DC, and neutrophils) and adaptive immune cells (e.g., CD4 + T cells, CD8 + T cells, and B cells) are the most critical components of tumor immunity. These cells participate in eliminating cancer cells by triggering inflammatory responses to inhibit tumor growth [62–64]. To investigate the impact of EZH2 expression on the immune infiltration landscape of pan-cancer, the link between EZH2 expression and infiltration levels of six tumor-infiltrating immune cells (CD4 + T cells, CD8 + T cells, B cells, macrophages, DC, and neutrophils) was assessed based on the TIMER algorithm. The results indicated that EZH2 expression was inversely or not correlated with these immune cells in most tumors, except for THYM (Figure 9A). The scatter plots of EZH2 expression and these immune infiltration cells in GBM, LUSC, TGCT, and sarcoma (SARC) further confirmed that EZH2 might suppress antitumor immunity by reducing the antitumor immune cell infiltration, leading to unfavorable prognosis in patients with various cancers (Figure 9C).

Furthermore, the correlation of EZH2 expression levels with infiltration of three immunosuppressive cells (MDSC, CAFs, and Tregs) was also examined. Positive correlations were observed between EZH2 expression and the infiltration of Tregs in BLCA, BRCA, BRCA-Basal, BRCA-LumA, CESC, GBM, HNSC, HNSC-HPV+, HNSC-HPV−, KIRP, LGG, LIHC, OV, PRAD, SKCM, STAD, and THCA; the infiltration of CAFs in ACC, CESC, KICH, LGG, LIHC, LUAD, MESO, PCPG, PRAD, SKCM, and THCA; the infiltration of MDSC in ACC, BLCA, BRCA, BRCA-Basal, BRCA-LumA, BRCA-LumB, CHOL, COAD, esophageal carcinoma (ESCA), GBM, HNSC-HPV−, KICH, KIRP, LGG, LIHC, LUAD, LUSC, MESO, OV, PAAD, PCPG, PRAD, READ, SARC, SKCM, STAD, TGCT, THYM, UCEC, UCS, and UVM (Figure 9B). It was also found that expression of EZH2 exhibited a strong positive correlation with dysfunctional T-cell phenotypes of KIRC, LUAD, BRCA-LumA, HNSC, myeloma, and colorectal cancer (Figure 9D).

In summary, EZH2 might facilitate tumor immune evasion by decreasing T-cell infiltration, but also promote T-cell dysfunction, ultimately resulting in the short survival time of cancer patients.

Immunotherapy is a new and booming treatment method that reactivates and restores the antitumor immune responses in the TME [65]. Given the significance of immune checkpoint molecule expression in predicting immunotherapy response, the correlation between EZH2 and expression levels of seven common immune checkpoint molecules (i.e., BTLA, CD276, LAG3, PD-1, PD-L1, PD-L2, and CTLA4) in TCGA pan-cancer cohort was tested to explore the potential value of EZH2 in immunotherapy [66]. The expression of EZH2 was positively correlated with the expression of immune checkpoint molecules in ACC, BLCA, CESC, CHOL, GBM, HNSC, KIRP, LAML, LGG, LGG, LUAD, LUSC, MESO, OV, PAAD, PCPG, READ, STAD, TGCT, THYM, UCEC, and UVM, LIHC, PRAD, especially in KIRC and THCA (Figure 10C). Besides that, EZH2 also exhibited a positive or negative correlation with the MHC and immunostimulatory genes in pan-cancer (Supplementary Figure S6). Thus, EZH2 might participate in the regulation of immune checkpoint molecules, affecting the pharmaceutical efficacy of the immune checkpoint inhibitor (ICI).

Frequent and constant mutations in tumor cells may contribute to clinical resistance against antitumor immunotherapy, leading to poor clinical outcomes in patients with cancers [67]. TMB, MSI, MMR, and neoantigens are associated with antitumor immune response in the TME, significantly affecting the response to immunotherapy [50,51,67,68]. Moreover, the evaluation results revealed that EZH2 expression was positively correlated with TMB in ACC, UCEC, THCA, STAD, SKCM, SARC, READ, PRAD, PAAD, OV, MESO, LUSC, LUAD, LIHC, LGG, KIRC, KICH, HNSC, COAD, BRCA, and BLCA but negatively correlated with TMB in THYM (Figure 10A). The correlation between MSI and EZH2 expression was positive in LUSC, LUAD, HNSC, GBM, BRCA, BLCA, UCEC, STAD, SARC, PRAD, MESO and negative in DLBC (Figure 10B). The functional defects in the MMR system caused irreparable errors in DNA replication, leading to frequent somatic mutations. After assessing the correlation between EZH2 expression and five MMR genes (i.e., MLH1, MSH2, MSH6, EPCAM, and PMS2), it can be seen that EZH2 expression was positively correlated to the expression levels of MMR genes in most cancers (Figure 10D). The results also indicated a significant positive correlation between EZH2 expression and neoantigens in LUAD (r=0.269, P=0.00048), STAD (r=0.441, P=9.09 × 10⁻¹³), HNSC (r=0.187, P=0.00175), BRCA (r=0.221, P=1.56 × 10⁻⁹), LGG (r=0.267, P=0.000152), and PRAD (r=0.307, P=5.09 × 10⁻⁷) (Figure 10E).

As reported in the previous literature, EZH2 activity and DNA methylation are closely associated with the resistance to immunotherapy [69]. Spearman correlation coefficients between EZH2 and four methyltransferase genes (DNMT1, 2, 3, 4) were calculated to investigate the correlation between EZH2 and DNA methylation. Subsequent results suggested that EZH2 was positively correlated with DNA methylation in almost all cancers, especially in LGG (r=0.81, P=5 × 10⁻¹²⁵), LIHC (r=0.8, P=1.4 × 10⁻⁸⁵) (Figure 10F).

Chemotherapy and targeted therapy have remained the most commonly adopted therapeutic approaches for various malignancies. Therefore, to address the value of EZH2 in chemotherapy and targeted therapy, the most sensitive and resistant antitumor agents targeting EZH2 were selected to guide clinical therapy selection based on the GDSC database. PARP-9482, PARP-0108, and Talazoparib were the most resistant anticancer drugs targeting EZH2 and may not be suitable for patients with EZH2 overexpression (Figure 11B). EZH2 expression had a significant positive correlation with susceptibility of Afuresertib and Venetoclax, which may be applied to patients with high EZH2 expressions (Figure 11A). Subsequent analyses indicated that Afuresertib and Venetoclax had higher IC50 values in the EZH2-wild subset, which may be more applicable to patients in the EZH2-mutation subset (Figure 11C,F). Furthermore, the structural formula and 3D structure of Afuresertib (C₁₈H₁₇Cl₂FN₄OS) were obtained through the DrugBank website (Figure 11D,E), and the molecular construct of Venetoclax (C₄₅H₅₀ClN₇O₇S) was displayed in Figure 11G.

To identify and characterize the molecular mechanism of EZH2 in tumorigenesis and development, a total of 20 EZH2-interacting genes (including EZH2) were selected through the GPS-Prot website to build a PPI network (Figure 12A). Then, functional pathway enrichment analysis (containing GO functional annotation and KEGG pathway enrichment analyses) was performed for these 20 EZH2 expression-associated genes. The results suggested that these genes were enriched in negative regulation of gene expression and epigenetic, ESC/E(Z) complex, PcG protein complex, histone methyltransferase complex, chromatin organization, methyltransferase complex, chromosome organization, histone lysine methylation, covalent chromatin modification, and peptidyl–lysine methylation (Figure 12B). For KEGG, the top ten pathways enriched were microRNAs in cancer, cysteine and methionine metabolism, nucleotide excision repair, longevity regulating pathway-multiple species, amphetamine addiction, cell cycle, ubiquitin-mediated proteolysis, cellular senescence, alcoholism, transcriptional misregulation in cancer (Figure 12C). To investigate the downstream signaling pathways affected by EZH2, the major biological functions of EZH2 and its molecular partners were detected by the GeneMANIA tool, mainly concentrated in PcG protein complex, regulation of gene expression and epigenetic, histone methyltransferase complex, methyltransferase complex, histone lysine methylation, G0 to G1 transition, regulation of G0 to G1 transition (Figure 12D). All samples in the TCGA pan-cancer cohort were classified into low‐ or high‐expression subgroups for GSEA analysis according to the median expression value of EZH2 (Figure 12E–H). The enrichment results of KEGG revealed that EZH2 expression was closely related to homologous recombination, cell cycle, nucleotide excision repair, asthma, complement and coagulation cascades, primary bile acid biosynthesis, and arachidonic acid metabolism. For Hallmark, EZH2 expression is mainly related to the G2M checkpoint, E2F targets, mTORC1 signaling, p53 pathway, myogenesis, bile acid metabolism, and coagulation.

To better investigate the upstream regulatory mechanisms of EZH2, Starbase and miRNet databases were adopted to screen out miRNAs targeting EZH2. Thirty-four miRNAs were predicted through the Starbase database, and 74 miRNAs were determined through the miRNet database (Supplementary Tables S1–S2). Based on the overlapping results of two websites, hsa-let-7a-5p and hsa-let-7b-5p were eventually identified, which were defined as the most significant miRNA regulators targeting EZH2 (Figure 12I). At the same time, the complementary sequences of EZH2 and target miRNAs (i.e., hsa-let-7a-5p and hsa-let-7b-5p) were shown in Figure 12J. Furthermore, 53 lncRNAs targeting hsa-let-7a-5p and hsa-let-7b-5p were also selected by the miRNet database to create an EZH2 ceRNA network (Figure 12K).

Since gastrointestinal malignancies are the most prevalent tumors, verifying the differential expression of EZH2 in digestive system cancers is the focus of the rest of the present study. First, our qRT-PCR results demonstrated that EZH2 was differentially expressed at the transcriptomic level in four gastrointestinal tumors, including LIHC, STAD, COAD, PAAD (Figure 13A–D). Moreover, we also performed the IHC experiments to compare the protein levels of EZH2 in various digestive system tumors. Subsequent IHC-staining images revealed that EZH2 remains differentially expressed in the four most common aerodigestive tract malignancies (Figure 13E–H).

As a major epigenetic regulator of transcriptional processes, the biological functions of EZH2 in autophagy, cell lineage determination, cell proliferation, DNA damage repair, and apoptosis have been demonstrated [70–73]. Its significant effect on the pathophysiology of cancer, including carcinoma occurrence, development, metastatic dissemination, immunomodulatory, and anticancer drug resistance, has also received widespread attention [25,26,74–80]. Since there is no pan-cancer study on EZH2, research on this topic is of great importance to deepen the understanding of cancer development and develop a novel therapeutic target for improving patient prognosis.

In the present study, genomics characteristics of EZH2 were first explored and analyzed. The human EZH2 gene maps on the seventh chromosome at region q36.1, including 20 exons and encoding 746 amino acids protein. Afterward, two domains of the EZH2 protein were investigated, namely CXC and SET domains.

The SET domain is the most predominant structure for maintaining the histone methyltransferase activity of EZH2. It is also critical for the formation of the N-terminus of the CXC domain [9]. EZH2 is conserved in various species and has maintained its major constituent structure (WD repeat-binging protein EZH2 and SET domain) during evolution. In addition, EZH2 could encode five isoforms involving the histone-lysine N-methyltransferase EZH2 isoforms a–e, mainly located in the nucleoplasm.

Then, the expression levels of EZH2 in different normal and tumor tissues were investigated. In normal human tissues, the mRNA expression levels of EZH2 in testis tissue were highest, followed by the esophagus and skin. The top three human cancer types with the most significant EZH2 mRNA expression were lung cancer, colorectal cancer, and head and neck cancer. The expression levels of EZH2 protein were roughly consistent with that of the mRNA levels in both normal and tumor tissues. Because of the few normal tissues in the TCGA database, the data from GTEx normal samples were integrated to examine the different expressions of EZH2 between tumor and corresponding paracancerous tissues. In addition to the significant down-regulation of LAML, EZH2 expression was significantly elevated in various cancer types, such as BLCA, BRCA, CESC, COAD, DLBC, GBM, LGG, LIHC, LUSC, OV, READ, STAD, THYM, UCEC, and UCS, consistent with previous experimental studies [17,25,81–88]. Further analysis also revealed that the expression of EZH2 was closely related to the pathologic stage or metastasis of LIHC, PRAD, KICH, KIRC, and lung cancer, suggesting that EZH2 could serve as one oncogene in tumor progression and metastasis.

The development and progression of cancer are caused by an accumulation of genetic alterations and epigenetic mutations [89]. Growing evidence indicates that genetic mutations can alter epigenetic patterns. Genetic modifications also cause instability and mutagenesis in the genome, demonstrating that genetic and epigenetic alterations may contribute to cancer formation [90–92]. Thus, the molecular features of EZH2 genetic alterations, DNA methylation, and protein phosphorylation were comprehensively detected and explored based on the TCGA pan-cancer cohort. In the TCGA pan-cancer cohort of 10953 patients, EZH2 mutations were detected in 290 cancer samples (2.6%), dominated by missense mutations and amplifications. Missense mutations were also the primary type of EZH2 SNV, with G > A (19.90%) accounting for the highest percentage. Moreover, EZH2 expression exhibited a positive relationship with CNV in numerous tumors, partially explaining the aberrant overexpression of EZH2 in most cancers. DNA methylation is one of the most prevalent epigenetic modifications, which can promote or inhibit gene expression [93]. Promoter hypermethylation usually suppresses gene expression, and hypomethylation means up-regulation of gene expression [94]. In the present study, the methylation status of EZH2 was lower in BLCA, UCEC, LUAD, LUSC, PRAD, READ, TGCT, and THCA than in normal tissues, and expression of EZH2 was up-regulated in these cancers, indicating that high expression of EZH2 might be caused by DNA methylation. Proteomic (phosphorylation) plays an essential role in investigating the molecular mechanisms of tumorigenesis and cancer progression [95,96]. According to the PhosphoSitePlus website, T487 was the most predominant phosphorylation site for EZH2 protein. The phosphorylation mechanisms of EZH2 protein in eight cancer types (GBM, LIHC, LUAD, breast cancer, PAAD, HNSC, UCEC, COAD, and OV) were explored using the CPTAC dataset.

GBM, LIHC, LUAD, and breast cancer had higher phosphorylation degrees in the T487 locus of EZH2 protein compared with normal tissues, suggesting that protein phosphorylation of EZH2 at the T487 site might function as a promoter in the development and progression of these cancers.

Subsequently, the association between EZH2 expression and prognosis was also evaluated. The overexpression of EZH2 was a risk factor for OS in ACC, KIRC, LGG, LIHC, MESO, PCPG, and PRAD patients. For THYM, the overexpression of EZH2 was a protective factor. Meanwhile, high expression of EZH2 was associated with poor DFS in patients with ACC, BLCA, KICH, KIRP, LGG, LIHC, PRAD, and THCA, indicating that EZH2 could serve as a predictive biomarker for pan-cancer prognosis.

The TME and immune escape correlate with tumor prognosis and treatment efficacy [97]. In general, tumor immune evasion involves two mechanisms. T-cell exclusion is the first mechanism of tumor immune escape by inhibiting the infiltration of immune cells, largely dependent on a high infiltration of immunosuppressive cells (i.e., MDSC, CAFs, and Tregs) [98,99]. In most tumors, EZH2 expression was inversely or not correlated with six immune cell types (CD4 + T cells, CD8 + T cells, B cells, macrophages, DC, and neutrophils) was positively correlated with the infiltration of Tregs, MDSC, and CAFs in multiple cancers. Second, T-cell dysfunction also accelerates the immune escape of tumor cells, resulting in tumor development, invasion, metastasis, and drug resistance [100,101]. According to the TIDE database, EZH2 expression showed a significant positive correlation with dysfunctional T-cell phenotypes of KIRC, LUAD, BRCA, HNSC, myeloma, and colorectal cancer. The above results demonstrated that EZH2 might facilitate tumor immune evasion by these two mechanisms (T-cell exclusion and dysfunction), leading to unfavorable survival outcomes in cancer patients.

To date, monoclonal antibodies targeting PD-1 or its ligand PD-L1 checkpoint appear to be the most effective and successful immunotherapeutic strategies available [102,103]. However, a few patients could benefit from PD-1/PD-L1 inhibitors [104,105]. Thus, it is imperative to develop inhibitors targeting novel immune checkpoint molecules to enhance the efficacy of existing ICI therapies [106]. In the present study, the correlation between EZH2 and the expression levels of seven common immune checkpoint molecules (i.e., BTLA, CD276, LAG3, PD-1, PD-L1, PD-L2, and CTLA4) was analyzed and examined. Further results suggested that EZH2 expression positively correlated to immune checkpoint molecule expression in most tumor types, especially in KIRC and THCA. Besides, as reported in the previous studies, TMB, MSI, MMR status, neoantigen, and methyltransferase gene were closely related to immunotherapy-associated response [52,67,107]. Spearman correlation analysis demonstrated a strong correlation between EZH2 expression and TMB, neoantigen, and MSI in multiple cancer types. MMR status and methyltransferase gene exhibited a coexpression relationship with EZH2. In conclusion, EZH2 had enormous potential as a new immunotherapy target.

Based on functional pathway enrichment analysis, it was found that EZH2 participated in tumor development and progression by modulating multiple signaling pathways and functions, such as negative regulation of gene expression and epigenetics and cell cycle. To investigate the upstream regulatory mechanisms of EZH2, two miRNA prediction authoritative databases (Starbase and miRNet) were used to predict miRNAs, and hsa-let-7a-5p and hsa-let-7b-5p were eventually identified as the most important miRNA regulators targeting EZH2 for constructing an EZH2 ceRNA network.

In the current study, the oncogenic mechanism of EZH2 in various cancers was elucidated by investigating the structure, expression, genetic mutation, DNA methylation, protein phosphorylation, and prognostic significance of EZH2 in the TCGA pan-cancer samples. Furthermore, EZH2 might participate in tumor immune evasion by T-cell exclusion and dysfunction. It exhibited a significant positive correlation with several important immunotherapy response biomarkers.

The datasets presented in the present study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

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

This work was supported by (the National Natural Science Foundation of China) (grant number 81860428), (the Key R & D general project of Jiangxi Science and Technology Department) (grant number 20203BBGL73187).

Lianghua Luo: Conceptualization, Resources, Data curation, Software, Formal Analysis, Supervision, Validation, Writing—original draft, Project administration, Writing—review & editing. Zhonghao Wang: Formal Analysis, Validation, Investigation, Visualization, Writing—original draft. Tengcheng Hu: Conceptualization, Resources, Data curation, Visualization. Zongfeng Feng: Conceptualization, Data curation, Software, Formal Analysis, Validation. Qingwen Zeng: Resources, Data curation, Software, Formal Analysis, Visualization. Xufeng Shu: Conceptualization, Resources, Data curation, Software, Formal Analysis. Ahao Wu: Resources, Data curation, Software. Pan Huang: Data curation, Software, Formal Analysis. Yi Cao: Supervision, Investigation. Yi Tu: Supervision, Funding acquisition, Writing—review & editing. Zhengrong Li: Supervision, Funding acquisition, Writing—original draft, Project administration, Writing—review & editing.

All contributors give consent for unrestricted publication of this work.

The studies involving human participants were reviewed and approved by Institutional Ethics Committee of the First Affiliated Hospital of Nanchang University. The patients/participants provided their written informed consent to participate in the prsent study.

Thanks to TCGA, GEO, GTEx, and CPTAC databases for providing gastric cancer sample data. TCGA, GEO, GTEx, and CPTAC belong to public databases. The patients involved in the database have obtained ethical approval. Users can download relevant data for free for research and publish relevant articles. Our study is based on open-source data, so there are no ethical issues and other conflicts of interest.

3D

three-dimensional

ACC

adrenocortical carcinoma

BLCA

bladder urothelial carcinoma

BRCA

breast invasive carcinoma

CAF

cancer-associated fibroblast

CDS

coding DNA sequences

ceRNA

competing endogenous RNA

CESC

cervical squamous cell carcinoma and endocervical adenocarcinoma

CHOL

cholangio carcinoma

CNA

copy number alteration

CNV

copy number variation

COAD

colon adenocarcinoma

COSMIC

catalogue of somatic mutations in cancer

CPTAC

Clinical proteomic tumor analysis consortium

CXC

cysteine-rich

DAB

diaminobenzidine

DC

dendritic cell

DFS

disease-free survival

DLBC

lymphoid neoplasm diffuse large B-cell lymphoma

DSS

disease-specific survival

EMT

epithelial–mesenchymal transition

ESCA

esophageal carcinoma

EZH2

enhancer of zeste homolog 2

GAPDH

glyceraldehyde 3-phosphate dehydrogenase

GBM

glioblastoma multiforme

GBMLGG

glioblastoma multiforme lower-grade glioma

GC

gastric cancer

GDSC

genomics of drug sensitivity in cancer

GEO

gene expression ontology

GO

gene ontology

GSEA

gene set enrichment analysis

GTEx

genotype-tissue expression

H3K27me3

trimethylation of Lys-27 in histone 3

HNSC

head and neck squamous cell carcinoma

HPA

human protein atlas

IBS

illustrator of biological sequences

IC50

half-maximal inhibitory concentration

ICI

immune checkpoint inhibitor

IHC

immunohistochemistry

K-M

Kaplan–Meier

KEGG

Kyoto Encyclopedia of Genes and Genomes

KICH

kidney chromophobe

KIPAN

pan-kidney cohort

KIRC

kidney renal clear cell carcinoma

KIRP

kidney renal papillary cell carcinoma

LAML

acute myeloid leukemia

LGG

brain lower grade glioma

LIHC

liver hepatocellular carcinoma

lncRNAs

long noncoding RNAs

LUAD

lung adenocarcinoma

LUSC

lung squamous cell carcinoma

m1A

N1-methyladenosine

m5C

5-methylcytosine

m6A

N6-methyladenosine

MDSC

myeloid-derived suppressor cell

MESO

mesothelioma

MMR

mismatch repair

MSI

microsatellite instability

MSI-H

microsatellite instability-high

MSI-L

microsatellite instability-low

NCBI

National Center for Biotechnology Information

OS

overall survival

OV

ovarian serous cystadenocarcinoma

PAAD

pancreatic adenocarcinoma

PcGs

polycomb group genes

PCPG

pheochromocytoma and Paraganglioma

PFS

progression-free survival

PPI

protein–protein interaction

PRAD

prostate adenocarcinoma

PRC2

polycomb repressive complex 2

qRT-PCR

quantitative real-time polymerase chain reaction

READ

rectum adenocarcinoma

SARC

sarcoma

SNV

single nucleotide variant

STAD

stomach adenocarcinoma

TCGA

The Cancer Genome Atlas

TGCT

testicular germ cell tumor

THCA

thyroid carcinoma

THYM

thymoma

TIDE

tumor immune dysfunction and exclusion

TMB

tumor mutation burden

Treg

regulatory T cell

UCEC

uterine corpus endometrial carcinoma

UCS

uterine carcinosarcoma

UVM

uveal melanoma