Multiomics characteristics and immunotherapeutic potential of EZH2 in pan-cancer

Abstract Enhancer of zeste homolog 2 (EZH2) is a significant epigenetic regulator that plays a critical role in the development and progression of cancer. However, the multiomics features and immunological effects of EZH2 in pan-cancer remain unclear. Transcriptome and clinical raw data of pan-cancer samples were acquired from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, and subsequent data analyses were conducted by using R software (version 4.1.0). Furthermore, numerous bioinformatics analysis databases also reapplied to comprehensively explore and elucidate the oncogenic mechanism and therapeutic potential of EZH2 from pan-cancer insight. Finally, quantitative reverse transcription polymerase chain reaction and immunohistochemical assays were performed to verify the differential expression of EZH2 gene in various cancers at the mRNA and protein levels. EZH2 was widely expressed in multiple normal and tumor tissues, predominantly located in the nucleoplasm. Compared with matched normal tissues, EZH2 was aberrantly expressed in most cancers either at the mRNA or protein level, which might be caused by genetic mutations, DNA methylation, and protein phosphorylation. Additionally, EZH2 expression was correlated with clinical prognosis, and its up-regulation usually indicated poor survival outcomes in cancer patients. Subsequent analysis revealed that EZH2 could promote tumor immune evasion through T-cell dysfunction and T-cell exclusion. Furthermore, expression of EZH2 exhibited a strong correlation with several immunotherapy-associated responses (i.e., immune checkpoint molecules, tumor mutation burden (TMB), microsatellite instability (MSI), mismatch repair (MMR) status, and neoantigens), suggesting that EZH2 appeared to be a novel target for evaluating the therapeutic efficacy of immunotherapy.


Introduction
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][3][4][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

Gene mutation characterization
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

Prognostic analysis
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.

Immune infiltration landscape
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].

Correlation analysis of EZH2 and several predictive biomarkers of immunotherapy
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][49][50][51][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.

The sensitivity and resistance analyses of anticancer agents targeting EZH2
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/).

Functional pathway enrichment analysis and construction of the ceRNA regulatory network
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 (lncR-NAs) targeting EZH2 miRNAs were determined by the miRNet database to create a competing endogenous RNA (ceRNA) regulatory network.

Experimentally validating the differential expression of EZH2
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 (CGGGAGCTGGAGCTAT-GATG); 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.

Genomics characteristics of EZH2
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).

Expression analysis of EZH2
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.

Mutational signature of EZH2
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 and protein phosphorylation of EZH2
DNA methylation patterns have been reported to contribute to tumorigenesis [58][59][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).

Prognostic analysis of EZH2
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.

The battle between tumor immunity and immunosuppressive microenvironment
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][63][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 infiltra-  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.

Prediction of anticancer drugs targeting EZH2
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 18 H 17 Cl 2 FN 4 OS) were obtained through the DrugBank website ( Figure 11D,E), and the molecular construct of Venetoclax (C 45 H 50 ClN 7 O 7 S) was displayed in Figure 11G.

Functional pathway enrichment analysis and construction of the ceRNA network
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).

Confirming the differential expression of EZH2
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).

Discussion
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][71][72][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][75][76][77][78][79][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][82][83][84][85][86][87][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][91][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.
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.

Conclusions
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.

Data Availability
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.