- Open Access
Epigenomic dysregulation-mediated alterations of key biological pathways and tumor immune evasion are hallmarks of gingivo-buccal oral cancer
Clinical Epigenetics volume 11, Article number: 178 (2019)
Gingivo-buccal oral squamous cell carcinoma (OSCC-GB) is the most common cancer among men in India and is associated with high mortality. Although OSCC-GB is known to be quite different from tongue cancer in its genomic presentation and its clinical behavior, it is treated identically as tongue cancer. Predictive markers of prognosis and therapy that are specific to OSCC-GB are, therefore, required. Although genomic drivers of OSCC-GB have been identified by whole exome and whole genome sequencing, no epigenome-wide study has been conducted in OSCC-GB; our study has filled this gap, and has discovered and validated epigenomic hallmarks of gingivobuccal oral cancer.
We have carried out integrative analysis of epigenomic (n = 87) and transcriptomic (n = 72) profiles of paired tumor-normal tissues collected from OSCC-GB patients from India. Genome-wide DNA methylation assays and RNA-sequencing were performed on high-throughput platforms (Illumina) using a half-sample of randomly selected patients to discover significantly differentially methylated probes (DMPs), which were validated on the remaining half-sample of patients.
About 200 genes showed significant inverse correlation between promoter methylation and expression, of which the most significant genes included genes that act as transcription factors and genes associated with other cancer types. Novel findings of this study include identification of (a) potential immunosuppressive effect in OSCC-GB due to significant promoter hypomethylation driven upregulation of CD274 and CD80, (b) significant dysregulation by epigenetic modification of DNMT3B (upregulation) and TET1 (downregulation); and (c) known drugs that can reverse the direction of dysregulation of gene expression caused by promoter methylation.
In OSCC-GB patients, there are significant alterations in expression of key genes that (a) regulate normal cell division by maintenance of balanced DNA methylation and transcription process, (b) maintain normal physiological signaling (PPAR, B cell receptor) and metabolism (arachidonic acid) pathways, and (c) provide immune protection against antigens, including tumor cells. These findings indicate novel therapeutic targets, including immunotherapeutic, for treatment of OSCC-GB.
Cancer of oral cavity and lip has the highest incidence (http://cancerindia.org.in/globocan-2018-india-factsheet/) and the highest mortality among men in India . About 90% of the oral malignancies are squamous cell carcinomas (SCC) — one of the most common histologic forms of head and neck squamous cell carcinoma (HNSCC) . Oral squamous cell carcinoma of the gingivo-buccal region (OSCC-GB), comprising buccal mucosa, gingivo-buccal sulcus, lower gingiva, and retromolar trigone, is the major subtype in India [3, 4]. Chewing of betel-quid containing areca nut (areca catechu) with or without tobacco, and also use of tobacco in other smokeless forms are prominent risk factors . Delayed presentation (~ 60% patients report at T3 or T4 stages) due to lack of awareness and a high rate of locoregional recurrence are major impediments to efficient management of oral cancer in India . Currently, the treatment of OSCC-GB is surgical resection followed by adjuvant radiotherapy with or without chemotherapy . There has been no significant decrease in the incidence of oral cancer for many decades . Despite aggressive multi-modality treatment, advanced oral cancer has a high rate of locoregional recurrence and poor prognosis with an overall 5-year survival rate of ~ 60%, which reduces further if associated with metastasis . There is a felt need to identify better prognostic and predictive markers to complement clinico-pathological findings in OSCC-GB.
In addition to genetic alterations, significant epigenetic alterations have been found to be associated with many cancer types . Promoter hypermethylation, mainly in CpG islands, can lead to downregulation of tumor suppressor genes and also dysregulation of downstream signaling pathways; these indications can serve as predictive markers of cancer . Global epigenomic alterations during the initiation and progression of cancer are indeed hallmarks of cancer .
We have previously identified a set of driver genes for OSCC-GB, with significantly enhanced burden of somatic mutations and copy number alterations, some of which are unique compared to drivers of other subtypes of head and neck cancers [9, 10]. Unique findings in OSCC-GB have also been reported by transcriptome-wide studies [11, 12]. OSCC-GB, therefore, deserves to be investigated separately. Differential methylation profiles of squamous cell carcinoma of different anatomic sites of the head and neck region have also been described [13,14,15,16,17,18], but none pertaining to the gingivo-buccal region of oral cavity. Among cancers of the oral cavity, tongue cancer is the most prevalent in western countries. Although OSCC-GB, the predominant type of oral cancer in India, is known to be quite different from tongue cancer in its molecular genetic presentation and its clinical behavior, it is treated identically as tongue cancer. Therefore, in addition to gaining deeper basic insights on epigenomic and transciptomic alterations in OSCC-GB, it is necessary to obtain predictive markers of prognosis and therapy that are specific to OSCC-GB.
Tumor and adjacent histologically normal tissue samples were collected in RNAlater, with written informed consent, from 101 OSCC-GB patients (details in Table 1), recruited at ACTREC, RADCH, and CNCI. TNM staging of each tumor sample was done following the 7th edition of the American Joint Committee on Cancer (AJCC) . Because the quality of biospecimens collected from the 101 patients varied, methylation and transcriptomic analyses could be performed on samples from 87 and 72 patients, respectively.
DNA extraction and bisulfite treatment
DNA was isolated using DNeasy Blood and Tissue Kit (QIAGEN). The purity and concentration was estimated using NanoDrop 2000 (Thermo Fisher Scientific). Approximately 500 ng genomic DNA from each sample was used for sodium bisulfite conversion using the EZ DNA methylation Gold Kit (Zymo Research, USA). Genome-wide DNA methylation was assayed using iScan (Illumina), for paired tumor and adjacent normal samples of 25 patients (all belonging to the “discovery” subset) using the Infinium Human Methylation450 BeadChip and of 62 patients (18 belonging to the “discovery” and 44 to “validation” subsets) using the Infinium MethylationEPIC BeadChip; these chips interrogate 485577 and 865918 CpG sites, respectively, of which 452512 CpG sites are common. For all assays, we followed the manufacturer’s protocols.
Processing of DNA methylation data
Raw array data (.IDAT files) were analyzed using the R package “minfi” and also using Illumina GenomeStudio Methylation module. The CpG-specific methylation level (β value), for each sample, was calculated as a ratio of fluorescent signal intensity of the methylated (m) and total of signal intensities from the methylated (m) and unmethylated (u) alleles as:
We have used the recommendation of the BeadChip manufacturer (Illumina) and set α = 100.
To reduce false positive inference, we have ignored any CpG probe (a) for which detection p was > 0.01; (b) “NA”- masked value; (c) that mapped to multiple locations on the human reference genome (hg19 with decoy sequence) when aligned using Bowtie 2 (v. 188.8.131.52) with “end-to-end” alignment mode and allowing for maximum 2 mismatches; (d) overlapped with a repetitive element [repeat masker (v. 4.0.5) (http://www.repeatmasker.org/) from UCSC hg19]; (e) polymorphic (MAF > 0.01) SNPs (dbSNP build 150) located within 10 bp of the interrogated CpG site; (f) spanned known regions of small insertions and deletions (indels) in the human genome (UCSC hg19); or (g) was located on a sex chromosome. Both intra-array (Infinium I and Infinium II) and inter-array (between samples) normalizations were performed: (a) subset quantile within array normalization (SWAN) and (b) inter-array quantile normalization (using β values).
Identification of epigenome-wide differential DNA methylation
Of the 87 patients, approximately 50% (43 patients, randomly selected) were used as a discovery cohort and the remaining (44 patients) were used as a validation cohort. We computed β values of the 452512 probes common to the two types of beadchip used for the 43 patients. β values of the additional 380341 probes for 18 patients whose samples were assayed using the MethylationEPIC BeadChip were computed and considered as “EPIC only” discovery-set. Wilcoxon signed rank test, corrected for multiple testing by the Benjamini-Hochberg method, was performed for each CpG site to test equality of distributions of β values between tumor and adjacent normal samples of patients. The CpG sites with Benjamini-Hochberg corrected p < 0.05 (merged set) or 0.02 (EPIC only set) and average |Δβ| ≥ 0.2 (tumor vs adjacent normal) were considered as significantly differentially methylated sites (DMPs). The stringency of the criteria used by us to identify significant DMPs is higher than those popularly used [13, 18]; this was done to avoid false positive inferences. A CpG site was considered hypermethylated if average Δβ was ≥ 0.2 or hypomethylated if average Δβ was ≤ − 0.2. Probes that were significantly differentially methylated in tumors compared with adjacent normal samples were used for the validation. A DMP was considered to be validated if for that DMP, in the validation cohort, the Benjamini-Hochberg corrected p was < 0.05 and the average |Δβ| was ≥ 0.2 (tumor vs adjacent normal). A significantly differentially methylated region (DMR) was defined as a gene region with a single or multiple unidirectional DMPs. The methylation level of a DMR was quantified as average of β values of DMPs that mapped to it.
RNA extraction and library preparation
RNA sequencing was performed on tumor and paired adjacent normal samples from 72 OSCC-GB patients. From each patient, total RNA was extracted from tumor and normal tissue samples using AllPrep DNA/RNA Mini Kit (QIAGEN). The quality and concentration of isolated total RNA were checked using Agilent 2100 Bioanalyzer and NanoDrop 2000 (Thermo Fisher Scientific). The OD260/OD280 ratio was ≥ 2, and RNA Integrity Number (RIN) was ≥ 7.0 for all sequenced samples. Ribosomal RNA (rRNA) was removed from the RNA samples, using Ribo-Zero Magnetic Kit (epicentre). Sequencing libraries were prepared from rRNA-depleted samples using TrueSeq RNA Sample Preparation Kit (Illumina). Each triplex cDNA library pool was sequenced as 100-bp paired-end on HiSeq-2000 or HiSeq-2500 (Illumina). Protocols suggested by the manufacturers were used for all assays.
Identification of transcriptome-wide differential gene expression
The 100-bp paired-end reads were aligned to the human reference transcript (build hg19) and reference genome (hg19 with decoy sequence) using the Tuxedo suite (TopHat2 (v. 2.1.1) , Bowtie 2 (v. 184.108.40.206) ), and SAMtools (v. 0.1.19) , with default settings. To minimize false positive inferences, multi-mapped and non-concordant reads were removed using SAMtools . Duplicate reads were identified using MarkDuplicates from PICARD (https://github.com/broadinstitute/picard) software (v 2.17.0) and discarded. Mapped transcripts were assembled using Cufflinks (v. 2.2.1) , with default parameters.
The set of 72 patients were randomly split into two equal subsets; these subsets were separately used for discovery and validation. Normalized gene expression (FPKM) values were estimated for the samples in both discovery and validation cohorts (each cohort n = 36) using cuffnorm. In the discovery cohort (n = 36), differential gene expression analysis between tumor and normal pairs was performed using cuffdiff with default values. Genes discovered to be significantly differentially expressed were considered for validation. Such a gene was declared as validated if the paired sample t test to compare mean expression levels in tumor and normal was significant (p < 0.05, after Benjamini-Hochberg multiple-testing correction).
Integrative analysis of DNA methylation and gene expression
To investigate the correlation between DNA methylation and gene expression in OSCC-GB patients, we have considered only those genes for which differential methylation in the promoter and differential expression level were both statistically significant. Such a gene was further examined to ascertain whether there was an opposite relationship between mean values of promoter methylation and gene expression in tumors compared to normals (i.e., hypomethylation and overexpression, or hypermethylation and underexpression). A validated gene that did not satisfy this property was discarded from further analysis. For each accepted gene, the Spearman correlation coefficient (ρ) was then calculated between the methylation level of DMRs (β values) and expression (FPKM) values observed in the patients. Genes for which the correlation coefficients were significant (Benjamini-Hochberg multiple-testing corrected p < 0.05) were used in pathway enrichment analysis.
Pathway enrichment analysis
A gene found to be significantly dysregulated and differentially methylated in the promoter was functionally annotated using the ClueGO plugin of Cytoscape (v. 3.7.0) . Gene Ontology (GO) terms and KEGG pathways were used in this analysis. To obtain significantly enriched categories of the GO terms and KEGG pathways, right-sided hypergeometric tests were performed; p values were corrected for multiple-testing using the Benjamini-Hochberg method.
Relationship with drug-induced gene expression
To investigate whether expression levels of genes found epigenetically altered can be altered in the opposite direction by a known drug, we have mined the Comparative Toxicogenomics Database (http://ctdbase.org/) of the DSigDB conglomerate database (http://tanlab.ucdenver.edu/DSigDB/DSigDBv1.0/).
DNA methylation and expression patterns
Epigenome-wide DNA methylation, assayed on paired tumor and adjacent normal tissue samples from 43 OSCC-GB patients (discovery cohort) and validated in an independent set of 44 patients (validation cohort), identified 25321 significant DMPs, of which 20023 were validated, of which 11387 were present on both Infinium 450 K and EPIC BeadChip arrays, while the remaining 8636 DMPs were present on EPIC BeadChip array only. Of these 20023 validated DMPs, 11522 (~ 58%) were hypermethylated and 8501 (~ 42%) were hypomethylated (Fig. 1). There is considerable variation in the proportion of DMPs across chromosomes; this proportion was highest for chromosome 8 and lowest for chromosome 16. With the exception of chromosomes 8 and 21, on all other chromosomes, the proportion of hypermethylated DMPs was higher than those hypomethylated (Fig. 1a); particularly high proportions of hypermethylated DMPs were observed for chromosomes 19 (82%), 17 (75%), and 22 (74%). Distinctly different patterns were observed between CGI (CpG islands) and non-CGI regions of the genome. Nearly half of the DMPs (~ 54%; 10744 of 20023) were enriched in non-CGI, open sea region of the genome. While about two-thirds of the DMPs in non-CGI regions were hypomethylated, almost all (~ 99%) DMPs located in the CGIs (~ 27% of all DMPs) were hypermethylated (Fig. 1b, c).
About 23% (4653 of 20023) of all DMPs mapped to promoter regions [1500 nt upstream of transcription start site (TSS) and from TSS to first exon] of known protein-coding genes. DNA methylation profiling identified and validated 6156 gene-regions that map to 4861 protein-coding genes. (It may be noted that multiple gene-regions can map to the same gene.) Significant differential methylation was observed in the promoter regions of 2318 genes, of which 1606 were hypermethylated and 712 were hypomethylated. About 38% (= 2318 of 6156) differentially methylated gene-regions belong to the promoter region of genes (Fig. 1d, e).
We discovered 2207 protein-coding genes to be significantly differentially expressed, at least by two-fold (mean of |log2 (fold change)| > 1), in tumors compared with paired normal samples, of which 1734 (78.6%) genes [804 (~ 46%) upregulated and 930 (~ 54%) downregulated] were validated (corrected p < 0.05; mean of |log2 (fold change)| > 1). Of the 1734 significantly differentially expressed genes, the maximum number of genes were on chromosome 1 (184 genes), followed by chromosomes 2 (127), 19 (125), 3 (110), and 11 (109).
Significant alteration of gene expression by epigenetic alteration associated with OSCC-GB
We validated that 307 genes were both significantly differentially methylated in promoters and also significantly differentially expressed. Of these, 209 genes were found to have significant negative correlation between β and FPKM values (Additional file 1: Table S1). Of these 209 genes, epigenetic downregulation due to significant promoter methylation was observed in 156 genes. The 10 genes with strongest negative correlation included ZNF132, ZNF626, ZSCAN18, ZNF844, SH2D2A, PHYHD1, IGF2BP2, ZNF829, ZNF880, and ZNF229 (Additional file 1: Table S1). Except SH2D2A and IGF2BP2, all genes were epigenetically downregulated.
Unsupervised hierarchical cluster analysis by Ward’s method using the Δβ values of 209 significantly differentially methylated genes identified two major clusters of patients (Fig. 2), comprising 15 and 29 patients, respectively. We note that the proportion of patients (55%) in the larger cluster who belongs to higher stages (T3 and T4) is greater compared to patients in the smaller cluster (33%). Further, about a third of the patients in the larger cluster presented with lymph node metastasis (N2); none of the patients in the smaller cluster presented with metastasis (Fig. 2).
Drug-induced reversal of epigenetic dysregulation
Of the 209 genes that we identified to be epigenetically dysregulated, the expression levels of 148 genes (Additional file 2: Table S2) are altered by known drugs in the direction opposite to that of epigenetic alteration; the total number of such expression-altering drugs is ~ 350. Of all significantly epigenetically downregulated (upregulated) genes in OSCC-GB, known drugs could upregulate (downregulate) 105 (43) genes.
Enriched pathways in OSCC-GB patients
Pathway enrichment analyses using 209 significantly dysregulated genes associated with epigenetic modifications in their promoter regions identified significant (corrected p < 0.05) enrichment of five KEGG pathways (Table 2 and Fig. 3) and twenty-two GO terms (Additional file 3: Table S3). Significantly enriched KEGG pathways included PPAR signaling (earlier found with significant gene expression alterations in OSCC-GB  and OSCC-tongue ), arachidonic acid metabolism (earlier implicated in OSCC using transcriptomic data  and in OSCC-GB on post-treatment disease-free survival length using somatic mutation data ), and B cell receptor signaling pathway (commonly implicated in chronic lymphocytic leukemia , but also in OSCC ). The enriched GO terms included negative regulation of G0 to G1 transition, process related to apoptosis such as regulation of leukocyte apoptotic process or negative regulation of dendritic cell apoptotic process, chemokine-mediated signaling pathway, and many other immune-related processes, including natural killer cell-mediated cytotoxicity, type I interferon signaling pathway. (Additional file 3: Table S3). The negative regulation of the cell cycle during G0 and G1 phases was shown in a previous study on OSCC . The chemokine-mediated signaling pathway was earlier found altered in lung cancer .
Comparison with TCGA data
The genome-wide DNA methylation profile observed in the OSCC-GB patients was validated using publicly available TCGA HNSCC data. In TCGA, all sub sites of HNSCC were included. There is considerable heterogeneity in regional presentation within the oral cavity . We have, therefore, considered data pertaining to the 31 OSCC-GB patients in the TCGA dataset (raw .IDAT files from the GDC Legacy Archive). In TCGA, OSCC-GB data were generated using the HM-450 K array while samples of our validation cohort patients were assayed using the HM-EPIC array platform. We, therefore, considered only those DMPs (CpG mapped to the promoter region of genes) that were common (80.4% of all promoter DMPs) to both arrays. Spearman’s correlation between β values of the 3741 common promoter DMPs of TCGA (n = 31) and our data (n = 44) of OSCC-GB patients was very high [ρ = 0.8574; p < 2.2 × 10−16] (Additional file 4: Figure S1). Patterns of clustering of patients based on β values of the common promoter DMPs were, therefore, largely similar (Additional file 4: Figure S2), thereby providing greater confidence to the inferences drawn from our data.
Gingivo-buccal oral squamous cell carcinoma, an anatomical and clinical subtype of head and neck squamous cell carcinoma (HNSCC) is the leading cancer among men in India. Epigenomic reprogramming, more specifically epigenomic gene repression, is a significant hallmark of oral cancer. Epigenome-wide DNA methylation profiling of 87 OSCC-GB patients (43 patients used as a discovery cohort and the remaining 44 as the validation cohort), using DNA methylation microarrays, identified 20023 significantly differentially methylated tumor-specific DMPs that mapped to 4861 protein-coding genes. Of these genes, 1606 were hyper- and 712 hypo-methylated in the promoter regions. Gene expression was significantly negatively correlated with DNA methylation in the promoter regions of 209 genes, of which 156 (~ 75%) were found epigenetically downregulated.
The ten genes with that showed high significant negative correlation between promoter methylation level and gene expression included mostly (7 of 10) zinc finger protein genes (ZNF132, ZNF626, ZSCAN18, ZNF844, ZNF829, ZNF880, and ZNF229) located on chromosome 19 (19q13; ZNF626 on 19p12). These genes encode C2H2-type of zinc finger proteins. Except ZNF880, all six zinc finger proteins are Krüppel-type family of transcription factors. A cluster of Krüppel-type zinc finger protein genes (including ZNF132 and ZSCAN18, earlier known as ZNF447) on chromosome 19q13 was earlier found significantly epigenetically downregulated in oropharyngeal squamous cell carcinoma . The remaining three genes are SH2D2A, PHYHD1, and IGF2BP2, of which IGF2BP2 is noteworthy. IGF2BP2 encodes an insulin-like growth factor-2 mRNA-binding protein and belongs to a conserved family of RNA-binding, oncofetal proteins. This protein modulates cell polarization, adhesion, and migration in tumor-derived cells and is strongly associated with cancer metastasis and the expression of oncogenic factors that include KRAS, MYC, and MDR1 . Upregulation of IGF2BP2 is associated with HNSCC . ZNF829 hypermethylation is associated with colorectal cancer . Overall, methylation of transcription factor and RNA binding genes and consequent alteration of expression seem important correlates of OSCC-GB.
Other significantly epigenetically downregulated genes are OSR1, SELENBP1, TGFBR3, and ZBTB16 [36,37,38,39]. Among significantly downregulated genes with promoter hypermethylation, CDON, ID4, and CPEB1 were associated with neuroblastoma, breast cancers, and hepatocellular carcinoma, respectively [40,41,42]. Upregulation with significant promoter hypomethylation of CD274, earlier known as PD-L1, was observed to be associated with OSCC-GB. Immunoexpression of the protein encoded by PD-L1 was found significantly increased in OSCC patients compared with normal or oral leukoplakia (oral pre-cancer) . Upregulation of another immune checkpoint gene CD80, observed in oral cancer , was found to be driven by significant promoter hypomethylation in OSCC-GB patients. CTLA4 that encodes an immune checkpoint receptor, was found to be transcribed in OSCC-GB tumour tissue by five- to six-fold (adjusted p < 0.05) higher than normal tissue (six-fold higher in discovery cohort; five-fold higher in validation cohort). However, this upregulation could not be ascribed to promoter hypomethylation of CTLA4 in our patients, contrary to the finding of a previous study on OSCC . Our observation of significant epigenetic upregulation of TRPM2 and LCK is consistent with earlier reports on oral cancer and oral lichen planus (oral pre-cancer) [44, 45]. We have also observed epigenetic upregulations of SULF1 and SEMA3C earlier found upregulated in gastric cancer  and breast cancer , respectively. However, our finding of upregulation of SEMA3C is not consistent with an earlier report on oral cancer . Our observations that DNMT3B, a gene responsible for de novo methylation process, is significantly upregulated with promoter hypomethylation and TET1, involved in DNA methylation process and gene activation, and is downregulated with promoter hypermethylation are worthy of further consideration.
We have also replicated some earlier observations on oral cancer and oral pre-cancer and in other cancers. These include epigenetic downregulation of GAS7, OSR1, SELENBP1, TGFBR3, and ZBTB16, earlier found associated with OSCC/HNSCC; epigenetic upregulation of TRPM2 and LCK (a proto-oncogene), earlier observed overexpressed in patients with OSCC or oral lichen planus (oral pre-cancer); altered expression of epigenetically downregulated (CDON, ID4, ZSCAN18, CPEB1, and NUPR1) or epigenetically upregulated (SULF1 and SEMA3C) genes earlier associated with several cancer types.
An important result of our study is that the immune checkpoint gene CD274 (earlier known as PD-L1) is significantly upregulated due to significant promoter hypomethylation. The protein encoded by PD-L1, a ligand of PD-1 receptor expressed on the surface of activated T cells, is expressed on the surface of antigen-presenting cells (APCs) which mostly include dendritic cells or macrophages. The induction of PD-L1 expression on the surface of tumor cells helps them evade immune response inhibition of the cytotoxic T cells and hence preventing apoptosis by the immunosuppressive properties of PD-L1/PD-1 interactions in the tumor microenvironment . Immunotherapy against immune checkpoints PD-L1 has shown promising results in treating metastatic HNSCC patients . Therefore, our finding prompts consideration of PD-1/PD-L1 immune checkpoint as a potential target for cancer therapy in OSCC-GB patients. The programmed death ligand 2 or PD-L2 (now known as PDCD1LG2) along with PD-L1 interacts with PD-1 to inhibit T cell activation and to downregulate the expression of anti-apoptotic molecules and the production of pro-inflammatory cytokines . We have observed significant upregulation of PDCD1LG2 in tumor samples. However, differential expression of this gene is driven by mechanisms other than differential methylation in OSCC-GB patients. In a recent study, it was shown that cis-PD-L1/CD80 interactions can disrupt PD-L1/PD-1 binding . This study  also observed that the binding was not affected by the simultaneous overexpression of CD80 and CD86. Similar to the previous study on OSCC patients , we have observed significant upregulation of CD80—driven by promoter hypomethylation—and CD86 in OSCC-GB patients. In addition, we have identified significant upregulation of CTLA4, and no significant alteration of expression or promoter DNA methylation of CD28 (which promotes the T cell activation). Therefore, in OSCC-GB patients, ligand-receptor interactions, such as, PD-L1/PD-1 or PD-L2/PD-1 and CTLA4/CD80 or CTLA4/CD86 may have crucial role in suppressing the effector T cell responses and hence immune evasion. These findings provide cues to further investigation that may be indicative of deployment of immune therapy in OSCC-GB.
We have noted upregulation of DNMT3B in OSCC-GB tumor samples driven by promoter hypomethylation. The product of DNMT3B acts as de novo DNA methyltransferase by transferring a methyl group from S-adenosylmethionine to the C-5 position of cytosine to form 5-methylcytosine (5mC), in the absence of a template . Overexpression of DNMT3B often correlates with the epigenetic inactivation of tumor suppressor genes leading to tumor formation , including oral tumor . On the contrary, the protein encoded by TET1 acts as demethylase which plays a key role in DNA demethylation by catalyzing the conversion of 5mC into 5-hydroxymethylcytosine (5hmC). It also mediates subsequent conversion of 5hmC into 5-formylcytosine (5fC) and 5-carboxylcytosine (5caC), which are further converted to unmethylated cytosine through the base excision repair pathway . Loss of TET1 expression through promoter CpG methylation frequently occurs in tumor cells, which results in tumor pathogenesis via inactivation of tumor suppressor genes . Epigenetic downregulation of TET1, observed in OSCC-GB patients, was earlier found associated with head and neck cancer . Therefore, our finding that epigenetic modifications upregulate DNMT3B and downregulate TET1, that likely plays an important role in OSCC-GB tumorigenesis, is novel.
We have identified that expression levels of about 71% (148 of 209) of the genes that are epigenetically dysregulated in OSCC-GB could be reversed by known drugs (Additional file 2: Table S2). The genes whose expression levels were altered by drugs included CD274, CD80, DNMT3B, TET1, PPARG, and PIK3CD. The drugs that induced reversal of expressions of these genes include dasatinib, calcitriol, tamoxifen, and aspirin. It is interesting that OSCC-GB caused by epigenetic alterations in important cancer genes is potentially actionable by known drugs.
Pathway analysis based on 209 significantly differentially genes identified PPAR signalling, arachidonic acid metabolism, and B cell receptor signaling pathways. The peroxisome proliferator-activated family of receptors (PPARs) are ligand-activated transcription factors. Epigenetic downregulation of PPARG (from the PPAR signaling pathway) and its ligands have earlier been shown to play a role in tumorigenesis. We have observed that PPARG is dysregulated in OSCC-GB patients and earlier also reported as dysregulated in HNSCC patients [57,58,59]. As a chemo-preventive agent, PPARG ligands act as targets for several types of cancer . Synthetic PPARG ligands have been shown to decrease incidence of carcinogen-induced tongue tumors in a dose-dependent manner . DNMT inhibitors (DNMTIs) are a promising class of anti-cancer therapies that renew transcription of a previously silenced gene . As a matter of fact, DNMTIs are being used to reverse epigenetic downregulation of PPARG, instead of PPARG ligands that show cytotoxicity. We have earlier shown  that survival benefits accrue to OSCC-GB patients who carry mutations in genes of the arachidonic acid metabolism (AAM) pathway. Our present finding of downregulation of oxidative stress-related genes EPHX2 and GPX3, belonging to the AAM pathway, in OSCC-GB patients, is therefore notable; similar findings were also reported in cancer of the oral cavity . Some genes that we found to be epigenetically upregulated in OSCC-GB patients, such as PIK3CD and IFITM1, were earlier reported to be upregulated in thyroid cancer  and HNSCC . IFITM1 plays an important role in invasion at early stages of HNSCC .
In OSCC-GB patients, we have observed a strong positive correlation of methylation levels of promoter DMPs between our study and the TCGA HNSCC study. Clusters of OSCC-GB patients of this study and patients included in TCGA, based on promoter DMPs that were common between the two studies, were similar. This is indicative of the robustness of the inferences of this study.
In sum, we have identified epigenetic alterations of key genes that regulate the feedback process of DNA methylation for the maintenance of normal cell division. Epigenetic alteration of transcription factor genes and RNA-binding genes and consequently expression dysregulation were strongly associated with OSCC-GB. Epigenetic activation of immunosuppressive PD-L1/PD-1 and CD80/CTLA4 interactions, which leads to inhibition of the cytotoxic T cells-mediated apoptotic process, was observed in gingivo-buccal oral cancer patients. Signaling (PPAR and B cell receptor) and the arachidonic acid metabolism pathways were found enriched in OSCC-GB. Our results provide novel therapeutic targets, including immunotherapy, for treatment of gingivo-buccal oral squamous cell carcinoma.
Availability of data and materials
Raw IDAT files of 174 samples generated using Illumina Infinium methylation array were deposited under EGAS00001003896 EGA study ID and aligned bam files for transcriptome data of 144 samples were deposited under EGAS00001003893 EGA study ID. Biospecimens may be shared on request, if available. The requester has to obtain prior approval from the Government of India to obtain biospecimens.
Differentially methylated probe
Differentially methylated region
Fragments per kilobase of transcript per million mapped reads
Head and neck squamous cell carcinoma
Kyoto Encyclopedia of Genes and Genomes
Minor allele frequency
Oral squamous cell carcinoma
Gingivo-buccal oral squamous cell carcinoma
Peroxisome proliferator-activated receptor
RNA Integrity Number
Squamous cell carcinomas
Subset quantile within array normalization
The Cancer Genome Atlas
Transcription start site
Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68:394–424.
Zanoni DK, Montero PH, Migliacci JC, Shah JP, Wong RJ, Ganly I, et al. Survival outcomes after treatment of cancer of the oral cavity (1985-2015). Oral Oncol. 2019;90:115–21.
Muttagi SS, Patil BR, Godhi AS, Arora DK, Hallikerimath SR, Kale AD. Clinico-pathological factors affecting lymph node yield in Indian patients with locally advanced squamous cell carcinoma of mandibular Gingivo-Buccal sulcus. Indian J Cancer. 2016;53:239–43.
Pathak KA, Gupta S, Talole S, Khanna V, Chaturvedi P, Deshpande MS, et al. Advanced squamous cell carcinoma of lower gingivobuccal complex: patterns of spread and failure. Head Neck. 2005;27:597–602.
More Y, D'Cruz AK. Oral cancer: review of current management strategies. Natl Med J India. 2013;26:152–8.
Klutstein M, Nejman D, Greenfield R, Cedar H. DNA methylation in cancer and aging. Cancer Res. 2016;76:3446–50.
Koch A, Joosten SC, Feng Z, de Ruijter TC, Draht MX, Melotte V, et al. Analysis of DNA methylation in cancer: location revisited. Nat Rev Clin Oncol. 2018;15:459–66.
Jones PA, Baylin SB. The epigenomics of cancer. Cell. 2007;128:683–92.
India Project Team of the International Cancer Genome Consortium. Mutational landscape of gingivo-buccal oral squamous cell carcinoma reveals new recurrently-mutated genes and molecular subgroups. Nat Commun. 2013;4:2873.
Bhosale PG, Cristea S, Ambatipudi S, Desai RS, Kumar R, Patil A, et al. Chromosomal alterations and gene expression changes associated with the progression of leukoplakia to advanced gingivobuccal cancer. Transl Oncol. 2017;10:396–409.
Prasad K, Rao R, Augustine D, Sowmya SV, Haragannavar V, Sagar P, et al. Pathway based prognostic gene expression profile of buccal and gingivo-buccal oral squamous cell carcinoma in smokeless tobacco chewers. Head Neck. 2019;41:388–97.
Singh R, De Sarkar N, Sarkar S, Roy R, Chattopadhyay E, Ray A, et al. Analysis of the whole transcriptome from gingivo-buccal squamous cell carcinoma reveals deregulated immune landscape and suggests targets for immunotherapy. PLoS One. 2017;12:e0183606.
Cancer Genome Atlas Network. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature. 2015;517:576–82.
Gaykalova DA, Zizkova V, Guo T, Tiscareno I, Wei Y, Vatapalli R, et al. Integrative computational analysis of transcriptional and epigenetic alterations implicates DTX1 as a putative tumor suppressor gene in HNSCC. Oncotarget. 2017;8:15349–63.
Krishnan NM, Dhas K, Nair J, Palve V, Bagwan J, Siddappa G, et al. A minimal DNA methylation signature in oral tongue squamous cell carcinoma links altered methylation with tumor attributes. Mol Cancer Res. 2016;14:805–19.
Pickering CR, Zhang J, Yoo SY, Bengtsson L, Moorthy S, Neskey DM, et al. Integrative genomic characterization of oral squamous cell carcinoma identifies frequent somatic drivers. Cancer Discov. 2013;3:770–81.
Chen Y, Hei N, Zhao J, Peng S, Yang K, Chen H, et al. A two-CpG-based prognostic signature for oral squamous cell carcinoma overall survival. J Cell Biochem. 2018. https://doi.org/10.1002/jcb.28182.
Basu B, Chakraborty J, Chandra A, Katarkar A, Baldevbhai JRK, Dhar Chowdhury D, et al. Genome-wide DNA methylation profile identified a unique set of differentially methylated immune genes in oral squamous cell carcinoma patients in India. Clin Epigenetics. 2017;9:13.
Sobin LH, Gospodarowicz MK, Wittekind C. TNM classification of malignant tumors. 7. Oxford: Wiley-Blackwell; 2010.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14:R36.
Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7:562–78.
Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics. 2009;25:1091–3.
Krishnan N, Gupta S, Palve V, Varghese L, Pattnaik S, Jain P, et al. Integrated analysis of oral tongue squamous cell carcinoma identifies key variants and pathways linked to risk habits, HPV, clinical parameters and tumor recurrence. F1000Res. 2015;4:1215.
Li S, Chen X, Liu X, Yu Y, Pan H, Haak R, et al. Complex integrated analysis of lncRNAs-miRNAs-mRNAs in oral squamous cell carcinoma. Oral Oncol. 2017;73:1–9.
Biswas NK, Das S, Maitra A, Sarin R, Majumder PP. Somatic mutations in arachidonic acid metabolism pathway genes enhance oral cancer post-treatment disease-free survival. Nat Commun. 2014;5:5835.
Ten Hacken E, Gounari M, Ghia P, Burger JA. The importance of B cell receptor isotypes and stereotypes in chronic lymphocytic leukemia. Leukemia. 2019;33:287–98.
Zhao E, Xu J, Yin X, Sun Y, Shi J, Li X. Detection of deregulated pathways to lymphatic metastasis in oral squamous cell carcinoma. Pathol Oncol Res. 2009;15:217–23.
Kudo Y, Kitajima S, Ogawa I, Miyauchi M, Takata T. Down-regulation of Cdk inhibitor p27 in oral squamous cell carcinoma. Oral Oncol. 2005;41:105–16.
Su H, Lin Z, Peng W, Hu Z. Identification of potential biomarkers of lung adenocarcinoma brain metastases via microarray analysis of cDNA expression profiles. Oncol Lett. 2019;17:2228–36.
Lleras RA, Adrien LR, Smith RV, Brown B, Jivraj N, Keller C, et al. Hypermethylation of a cluster of Krüppel-type zinc finger protein genes on chromosome 19q13 in oropharyngeal squamous cell carcinoma. Am J Pathol. 2011;178:1965–74.
Bell JL, Wächter K, Mühleck B, Pazaitis N, Köhn M, Lederer M, et al. Insulin-like growth factor 2 mRNA-binding proteins (iGF2BPs): post-transcriptional drivers of cancer progression? Cell Mol Life Sci. 2013;70:2657–75.
Davidson MA, Shanks EJ. 3q26-29 Amplification in head and neck squamous cell carcinoma: a review of established and prospective oncogenes. FEBS J. 2017;284:2705–31.
Naumov VA, Generozov EV, Zaharjevskaya NB, Matushkina DS, Larin AK, Chernyshov SV, et al. Genome-scale analysis of DNA methylation in colorectal cancer using Infinium HumanMethylation450 BeadChips. Epigenetics. 2013;8:921–34.
Chen W, Zou A, Zhang H, Fu X, Yao F, Yang A. Odd-skipped related transcription factor 1 (oSR1) suppresses tongue squamous cell carcinoma migration and invasion through inhibiting NF-κB pathway. Eur J Pharmacol. 2018;839:33–9.
Chen F, Chen C, Qu Y, Xiang H, Ai Q, Yang F, et al. Selenium-binding protein 1 in head and neck cancer is low-expression and associates with the prognosis of nasopharyngeal carcinoma. Medicine (baltimore). 2016;95:e4592.
Meng W, Xia Q, Wu L, Chen S, He X, Zhang L, et al. Downregulation of TGF-beta receptor types II and III in oral squamous cell carcinoma and oral carcinoma-associated fibroblasts. BMC Cancer. 2011;11:88.
Boldrup L, Gu X, Coates PJ, Norberg-Spaak L, Fahraeus R, Laurell G, et al. Gene expression changes in tumor free tongue tissue adjacent to tongue squamous cell carcinoma. Oncotarget. 2017;8:19389–402.
Gibert B, Delloye-Bourgeois C, Gattolliat CH, Meurette O, Le Guernevel S, Fombonne J, et al. Regulation by miR181 family of the dependence receptor CDON tumor suppressive activity in neuroblastoma. J Natl Cancer Inst. 2014;106.
Nasif D, Campoy E, Laurito S, Branham R, Urrutia G, Roqué M, et al. Epigenetic regulation of ID4 in breast cancer: tumor suppressor or oncogene? Clin Epigenetics. 2018;10:111.
Xu M, Fang S, Song J, Chen M, Zhang Q, Weng Q, et al. CPEB1 mediates hepatocellular carcinoma cancer stemness and chemoresistance. Cell Death Dis. 2018;9:957.
Stasikowska-Kanicka O, Wągrowska-Danilewicz M, Danilewicz M. CD8+ and CD163+ infiltrating cells and PD-L1 immunoexpression in oral leukoplakia and oral carcinoma. APMIS. 2018;126:732–8.
Zhao LY, Xu WL, Xu ZQ, Qi C, Li Y, Cheng J, et al. The overexpressed functional transient receptor potential channel TRPM2 in oral squamous cell carcinoma. Sci Rep. 2016;6:38471.
Oluwadara O, Giacomelli L, Christensen R, Kossan G, Avezova R, Chiappelli F. LCK, survivin and PI-3 K in the molecular biomarker profiling of oral lichen planus and oral squamous cell carcinoma. Bioinformation. 2009;4:249–57.
Junnila S, Kokkola A, Mizuguchi T, Hirata K, Karjalainen-Lindsberg ML, Puolakkainen P, et al. Gene expression analysis identifies over-expression of CXCL1, SPARC, SPP1, and SULF1 in gastric cancer. Genes Chromosomes Cancer. 2010;49:28–39.
Cole-Healy Z, Vergani P, Hunter K, Brown NJ, Reed MW, Staton CA. The relationship between semaphorin 3C and microvessel density in the progression of breast and oral neoplasia. Exp Mol Pathol. 2015;99:19–24.
Dermani FK, Samadi P, Rahmani G, Kohlan AK, Najafi R. PD-1/PD-L1 immune checkpoint: Potential target for cancer therapy. J Cell Physiol. 2019;234:1313–25.
Ghanizada M, Jakobsen KK, Grønhøj C, von Buchwald C. The effects of checkpoint inhibition on head and neck squamous cell carcinoma: a systematic review. Oral Oncol. 2019;90:67–73.
Sugiura D, Maruhashi T, Okazaki IM, Shimizu K, Maeda TK, Takemoto T, et al. Restriction of PD-1 function by cis-PD-L1/CD80 interactions is required for optimal T cell responses. Science. 2019;364(6440):558–66.
Baubec T, Colombo DF, Wirbelauer C, Schmidt J, Burger L, Krebs AR, et al. Genomic profiling of DNA methyltransferases reveals a role for DNMT3B in genic methylation. Nature. 2015;520:243–7.
Zhang Y, Charlton J, Karnik R, Beerman I, Smith ZD, Gu H, et al. Targets and genomic constraints of ectopic Dnmt3b expression. Elife. 2018;7.
Chen WC, Chen MF, Lin PY. Significance of DNMT3b in oral cancer. PLoS One. 2014;9:e89956.
Wu X, Zhang Y. TET-mediated active DNA demethylation: mechanism, function and beyond. Nat Rev Genet. 2017;18:517–34.
Li L, Li C, Mao H, Du Z, Chan WY, Murray P, et al. Epigenetic inactivation of the CpG demethylase TET1 as a DNA methylation feedback loop in human cancers. Sci Rep. 2016;6:26591.
Misawa K, Imai A, Mochizuki D, Mima M, Endo S, Misawa Y, et al. Association of TET3 epigenetic inactivation with head and neck cancer. Oncotarget. 2018;9:24480–93.
Elrod HA, Sun SY. PPARgamma and apoptosis in cancer. PPAR Res. 2008;2008:704165.
Burotto M, Szabo E. PPARγ in head and neck cancer prevention. Oral Oncol. 2014;50:924–9.
Yoshida K, Hirose Y, Tanaka T, Yamada Y, Kuno T, Kohno H, et al. Inhibitory effects of troglitazone, a peroxisome proliferator-activated receptor gamma ligand, in rat tongue carcinogenesis initiated with 4-nitroquinoline 1-oxide. Cancer Sci. 2003;94:365–71.
Mack GS. To selectivity and beyond. Nat Biotechnol. 2010;28:1259–66.
Pedro NF, Biselli JM, Maniglia JV, Santi-Neto D, Pavarino ÉC, Goloni-Bertollo EM, et al. Candidate biomarkers for oral squamous cell carcinoma: differential expression of oxidative stress-related genes. Asian Pac J Cancer Prev. 2018;19:1343–9.
Bu Q, You F, Pan G, Yuan Q, Cui T, Hao L, et al. MiR-125b inhibits anaplastic thyroid cancer cell migration and invasion by targeting PIK3CD. Biomed Pharmacother. 2017;88:443–8.
Hatano H, Kudo Y, Ogawa I, Tsunematsu T, Kikuchi A, Abiko Y, et al. IFN-induced transmembrane protein 1 promotes invasion at early stage of head and neck cancer progression. Clin Cancer Res. 2008;14:6097–105.
We are grateful to all participating members of Systems Medicine Cluster (SyMeC) and International Cancer Genome Consortium (ICGC) India project for their inputs and contributions to this study.
The present work was supported by grant from Department of Biotechnology (DBT), Govt. of India, through the SyMeC project (BT/Med-II/NIBMG/SyMeC/2014/Vol. II). Patient recruitment and sample collection in ACTREC, Mumbai, were supported by funds provided by DBT to the ICGC’s India project.
Ethics approval and consent to participate
The study was approved by the Institutional Ethics Committees of Advanced Centre for Treatment, Research and Education in Cancer (ACTREC), Mumbai, Dr. R. Ahmed Dental College & Hospital (RADCH), Kolkata, Chittaranjan National Cancer Institute (CNCI), Kolkata, National Institute of Biomedical Genomics (NIBMG), Kalyani, and Indian Statistical Institute (ISI), Kolkata. Prior written informed consent was obtained from all study participants.
Consent for publication
All authors have agreed to publish in Clinical Epigenetics.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. Relevant details of 209 genes with significant negative correlation between promoter methylation and gene expression.
. Information related to alteration of expression of 209 epigenetically dysregulated genes by one or more known drugs.
Significantly enriched GO terms in OSCC-GB patients based on 209 genes significantly differentially methylated in their promoter regions and related information.
. Scatter diagram showing the relationship of average β values of the DMPs found in the present study with those in the TCGA study, irrespective of whether these probes were also significantly differentially methylated in the TCGA study. Each point on the scatter diagram indicates for a DMP of the present study the average β value over the 44 OSCC-GB patients included in the validation subset and, for the TCGA study, averaged over the 31 OSCC-GB patients. Figure S2. Integrated unsupervised hierarchical clustering and heatmap using quantitative methylation level of 3741 DMPs mapped to gene promoter from (A) 44 OSCC-GB patients included in this study and (B) 31 OSCC-GB patients included in TCGA study. Corresponding clusters obtained from two sets of patients showed similar methylation pattern and phenotypic features. [White box in top panels of (A) and (B) indicates unavailability of the respective clinical information].
About this article
Cite this article
Das, D., Ghosh, S., Maitra, A. et al. Epigenomic dysregulation-mediated alterations of key biological pathways and tumor immune evasion are hallmarks of gingivo-buccal oral cancer. Clin Epigenet 11, 178 (2019). https://doi.org/10.1186/s13148-019-0782-2
- Gingivo-buccal oral squamous cell carcinoma
- Integrative analysis
- Immunotherapeutic marker