Genome-Wide DNA Methylation Confirms Oral Squamous Cell Carcinomas in Proliferative Verrucous Leukoplakia as a Distinct Oral Cancer Subtype: A Case Control Study

1. Introduction

Proliferative verrucous leukoplakia (PVL) is the oral potentially malignant disorder (OPMD) with the highest transformation rate into oral cancer [1,2,3]. Recently, in an update coordinated by the World Health Organization Collaborating Centre for Oral Cancer, this disease of unknown origin has been officially recognized for the first time as a distinct form of classical oral leukoplakias (cOL) characterized by a progressive, persistent, and irreversible clinical course with the presence of multiple oral leukoplakias that frequently become verrucous [1]. In early stages, there are no pathognomonic clinicopathological features that allow differentiate from frictional keratosis, homogeneous leukoplakias, or oral lichenoid diseases [4,5]. In advanced stages, lesions progress to multifocality with a wide spectrum of lesion types ranging from smooth or fissured white patches to exophytic, verrucous, erythematous, and erythroleukoplakia lesions [6,7,8]. Finally, during their evolution, it is observed that between 43.87% (95% CI: 31.93 – 56.13) [3] to 65.8% (95% CI: 55.3 – 76.2) [2] of patients develop at least one malignant transformation into oral cancer, although some authors still indicate higher malignancy rates [5,6,7].

PVL, unlike classic oral leukoplakia, predominantly affects women at advanced ages, without tobacco or alcohol habits [1,7,9], and primarily occurs in different locations such as the gums, palate, alveolar mucosa, and buccal mucosa [1,6,7]. Moreover, the malignant transformation rate is five times higher [3] with high rates of recurrence [10] and challenging management.

Head and neck carcinomas comprise a heterogeneous and aggressive group of malignant neoplasm located in different areas of the upper aerodigestive tract. These carcinomas have different risk factors, molecular characteristics, clinical behaviour, and treatment response [11,12,13]. Oral squamous cell carcinoma (OSCC) is one of the most aggressive types characterized by poor prognosis and high rates of metastasis and locoregional recurrences [14]. For years, it has been discussed and hypothesized that oral cavity cancers present subtypes with clinical, biological, and evolutionary differences, such as tongue carcinomas, carcinomas associated with the human papillomavirus, or carcinomas derived from oral potentially malignant disorders [1,15,16]. Therefore, it is reasonable to assume that oral squamous cell carcinomas in patients with proliferative verrucous leukoplakia (PVL-OSCC) may represent a differentiated phenotype of oral cavity squamous cell carcinoma and this phenotype could be explained by differential epigenetics.

Recent reviews have summarised the molecular alterations described in PVL and OSCC [17,18]. Weak evidence has been suggested for chromosomal instability such as DNA aneuploidy, loss of heterozygosity at locus 9p21, or specific expression of Mcm protein in PVL [18]. All these genetic alterations ultimately result in aberrant gene expression. Although some transcriptomic studies have revealed differential expression genes in PVL related to immune surveillance [19] or with upregulation of many cancers associated genes [20], their aetiology and pathogenesis have not been fully defined [18]. The complexity of PVL disease and the landscape of genetic alterations is insufficient to explain the pervasive gene expression changes and alterations to cellular function in cancer. A new focus lies on gene expression regulatory events such as epigenetic factors. DNA methylation is one of the most well-studied epigenetic mechanisms, and it plays an essential role in regulating gene expression [21,22]. Using high-definition microarray approach is possible to evaluate differential DNA methylation patterns, covering a representative portion of all human CpG sequences throughout the genome [22].

Previous studies from our group [23,24] as well as other group [25] have elucidated the implication of aberrant DNA methylation in the pathogenesis of PVL. These studies confirmed the role of epigenetic alterations and the potential of methylation markers in PVL without history of oral cancer compared to healthy controls and suggested novel OSCC prognosis biomarkers [23]. Recently, in an integrative analysis of transcriptomic and methylomic data in patients with PVL-OSCC revealed that 20% of the 133 differentially expressed genes were potentially regulated by DNA methylation [24], suggesting differences in their transformation pathways [24,25].

Although carcinomas in PVL patients present different clinical, histopathological, and evolutionary features compared to conventional carcinomas, differences at molecular level have hardly been studied. The main objectives were (1) to compare the methylome profile differences between patients with oral squamous cell carcinoma preceded or not by proliferative verrucous leukoplakia in an epigenome-wide association study (EWAS) using the high-coverage Illumina Infinium HumanDNAMethylation 850 BeadChip and (2) to identify potential prognostic biomarkers.

2. Materials and Methods

• Study design and Patients

We conducted an EWAS on 18 patients treated at the Stomatology and Maxillofacial Surgery Department of Hospital General Universitario de Valencia. Participants provided written inform consent, and study protocols and procedures were approved in compliance with the Helsinki Declaration by the Ethics Committee for Human Research of the Universitat de València (Ref. H1523722754549) and the Consortium of the General University Hospital of Valencia (30 May 2019).

This case – control study was conducted in two patient groups: Group I included eight patients with the development of the fist oral squamous cell carcinoma arising with PVL (OSCC – PVL) and Group II included ten patients with conventional oral squamous cell carcinoma without an antecedent of any oral potentially malignant disorder (cOSCC). PVL was diagnosed according to Villa et al. (2018) criteria [7]. Clinical and pathological TNM stage was classified according to the 8th edition of the American Joint Commission on Cancer, 2017 [26]. For each patient two representative biopsies including epithelial and connective tissue from the cancerous area were obtained between 2019 and 2021 for the histopathological and methylation studies. The latter were immediately frozen at -80ºC until the analysis.

• DNA extraction and Genome – Wide DNA Methylation

Total DNA from fresh frozen tissue was extracted using the column-based DNA extraction method (E.Z.N.A. DNA kit and DNeasy Blood & Tissue Kit; Qiagen, Hilden, Germany) according to the manufacturer’s recommendations. DNA concentration and quality control was quantified by Qubit 2.0 Fluorometer (Quant-iT PicoGreen dsDNA Assay, Life Technologies, CA, USA). DNA integrity was evaluated by electrophoresis in a 1.3% agarose gel.

For methylation profiling, we used the validated Infinium HumanMethylationEPIC BeadChip (850K) array (Illumina Inc., San Diego, CA, USA) which interrogates over 862,927 CpG sites across the genome and covers 99% of known genes and 95% of CpG islands [27]. Only samples providing 500ng of high-quality DNA were used for bisulfite conversion using the EZ-96 DNA Methylation kit (Zymo Research Corp., CA, USA) following the manufacturer’s instructions for Infinium assays. A total of 4 μl of the bisulphite converted DNA were processed following the Illumina HD Methylation Assay Protocol. The Bead Chips were scanned with an Illumina iScan System and intensity values (.idat files) were generated.

• Data normalisation

Raw data (IDATs) were normalized using the minfi R package (v 1.38) and functional normalization. CpG markers present on MethylationEPIC were classified based on their chromosome location, the Infinium chemistry used to interrogate the marker (Infinium I, Infinium II), and the feature category gene region as per UCSC annotation (TSS200, TSS1500, 5′UTR, 1st Exon, Body, 3′UTR). Additional criteria included the location of the marker relative to the CpG island (open sea, island, shore, shelf). Probe filtering was performed at three levels. Each beta value in the EPIC array was accompanied by a detection p-value that represents the confidence of a given beta value. CpGs with a high p-value p > 0,01 (1,620 CpGs) were removed to exclude probes with a low quality and high variability (28). Probes overlapping with single nucleotide polymorphisms were removed because it can alter methylation levels (2,932 CpGs), as well as probes associated with sex chromosomes (19,681 CpGs) to avoid potential bias due to gender. Therefore, after filtering, 842,179 CpGs were considered valid for the study. For each CpG site, a specific β value was obtained that ranged from 0 (no methylation) to 1 (complete methylation). The difference in mean β-values between the groups was indicated as Δβ.

• Statistical analysis and bioinformatics

To summarize the characteristics of the different groups, descriptive statistical tests were performed. Chi-square tests were used to compare proportions. Student t-tests and ANOVA tests were applied to compare crude means of continuous variables. For all tests, p < 0.05 was considered to indicate statistical significance.

Firstly, the unsupervised exploratory analysis was performed using a principal component analysis (PCA) and a heatmap of claustered observations and variables where green represents a higher methylation gain. Secondly, supervised differential methylation analysis was assessed using two different approaches. On the one hand, using a rank-based regression model for each CpG [29]. p-values for the regression coefficients were adjusted for multiple comparisons by using False Discovery Rate. Adjusted p-values lower than 0.05 were considered statistically significant. On the other hand, an elastic net penalized logistic regression model was adjusted to select the CpGs able to discriminate between groups [30]. The penalization factor for the elastic net was selected by taking the highest lambda at one standard error from the minimum (one-standard-error rule) from 500 repetitions of 10-fold cross-validation. Then, the median of the 500 lambda values was used as the final penalization factor. Elastic net alpha parameter (regulating the mix of L1 and L2 penalization) was set at 0.4. CpGs with non-zero coefficients after the penalization were selected and thus were considered as relevant for discriminating between both groups. All statistical analyses were performed using R (v 4.2.0), R packages glmnet (v 4.1-4) and Rfit (v 0.24.2).

3. Results

3.1. Patient Characteristics

The clinicopathological, demographic and evolutionary characteristics of patients are described in Table 1. The mean age and follow up was higher in PVL-OSCC patients. We only found significant differences in the clinical form of tumour presentation, being more variable in PVL-OSCC as they presented exophytic, erythroplastic or mixed lesions. Interestingly, PVL-OSCC showed a tendency towards a more favourable prognosis, characterized by less perineural infiltration, reduced depth of invasion, fewer occurrences of cervical lymph nodes metastasis and an earlier TNM stage (Table 1).

3.2. The Unsupervised Exploratory Analysis did Not Allow to Differentiate Between the Different Oral Cancer Groups

In this first approach, an exploratory principal component analysis didn‘t reveal clear differences between both groups (Figure 1A). As expected, no distinct clusters or outliers were observed among the different patient groups. Next, we performed an unsupervised heatmap using 5,000 randomly selected CpGs and performed hierarchical clustering analysis on both observations and variables (Figure 1B). The findings were consistent with the previous analysis, showing with no clear distinction between groups. This lack of distinct clusters can be explained by the fact that both groups consist of patients with oral squamous cell carcinomas.

3.3. Identification of Differentially Methylated CpG Sites

Two supervised differential methylation analysis were performed. The first one by adjusting a rank-based regression model for each CpG methylation value. After correcting for multiple comparisons using False Discovery Rate, 14 CpGs corresponding to 10 genes showed statistically significant differences between both groups with an adjusted p-value < 0,05 (Table 2). Additionally, the results were also represented in a heatmap with the values of these CpGs and the differences were visually confirmed using the 14 CpGs (Figure 2).

An additional second analysis was performed to define the existence of potential specific methylation differences using an elastic net-penalized logistic regression model. The 500 repetitions of the cross-validation procedure selected a lambda value of 0.95. This statistical analysis determined 9 CpGs which provide the highest discrimination power (Table 3 and Figure 3).

We also observed that 16 of the 23 CpGs (70%) were hypomethylated in PVL-OSCC group.

Table 3. Differentially methylated CpGs from the elastic net-penalized logistic regression model.

MethylationEPIC probe ID	Gene symbol	Cytoband	∆β	Relation to CpG island	Regulatory feature
cg24794734	AGL	1p21.2	-0.007	Shore	1st Exon, 5’UTR, TSS1500
cg11369761	-		0.12	Shore	-
cg15199741	FGD5	3p25.1	-0.09	Open sea	Body
cg09643539	ARL15	5q11.2	-0.14	Open sea	Body
cg09365002	DAXX	6p21.32	-0.3	Shore	Body
cg17251196			-0.17	Shore	Body
cg24104268	-		-0.12	Island	-
cg12003952	ZNF429	19p12	-0.07	Shore	TSS1500
cg15977137	WRB / GET1	21q22.2	-0.16	Shore	Body, 1st Exon, 5’UTR

Table 3. Differentially methylated CpGs from the elastic net-penalized logistic regression model.

MethylationEPIC probe ID	Gene symbol	Cytoband	∆β	Relation to CpG island	Regulatory feature
cg24794734	AGL	1p21.2	-0.007	Shore	1st Exon, 5’UTR, TSS1500
cg11369761	-		0.12	Shore	-
cg15199741	FGD5	3p25.1	-0.09	Open sea	Body
cg09643539	ARL15	5q11.2	-0.14	Open sea	Body
cg09365002	DAXX	6p21.32	-0.3	Shore	Body
cg17251196			-0.17	Shore	Body
cg24104268	-		-0.12	Island	-
cg12003952	ZNF429	19p12	-0.07	Shore	TSS1500
cg15977137	WRB / GET1	21q22.2	-0.16	Shore	Body, 1st Exon, 5’UTR

Figure 3. Hierarchical clustering and heatmap of the 9 CpGs predicted by the elastic net-penalized logistic regression model.

Figure 3. Hierarchical clustering and heatmap of the 9 CpGs predicted by the elastic net-penalized logistic regression model.

3.4. Identification of Potential DNA Methylation Biomarkers

Interestingly, upon crossing the results of the two supervised analysis, it was observed that there were two genes that were significant differentially methylated in both, one of them corresponding to the WRB gene and the other to the AGL gene.

Furthermore, we used data from a previous study published by our group [31] describing differential methylation patterns between PVL without antecedent or presence of oral cancer and homogeneous leukoplakias and healthy controls. An integrative analysis was conducted using the 163 differentially methylated CpGs described in PVL patients (none of which are part of the present study). ARL15 was the only gene differentially methylated in PVL patients with or without malignant transformation when compared to homogeneous leukoplakias, cOSCC and healthy controls.

4. Discussion

The aggressiveness of PVL, characterized by multifocality, high recurrence rates, as well as malignant transformation and a tendency to develop second primary tumours, suggests the presence of molecular alterations not yet characterized in PVL patients [10,18]. Here, we compare the genome wide DNA methylation signatures in fresh frozen tissues between conventional OSCC and PVL-OSCC in an attempt to understand whether epigenetic dysregulation might explain their clinical and evolutionary differences.

First of all, we observed that PVL-OSCC patients had a more favourable prognosis than cOSCC. This was evidenced by two published meta-analyses [32,33]), especially with regard to the mortality rate (21.29%. CI: 8.77 – 36.36)[32], recurrence rate (22%)[33] and the uncommon development of metastases in PVL[33]. In contrast, cOSCC exhibited higher mortality rates [34] and cervical lymph nodes metastasis (40%)[35]. We aim to comprehend if these distinct clinical phenotypes could be explained by differential methylome.

Studying epigenetic heterogeneity in cancer evolution is critical because the epigenetic dysregulation is a hallmark of cancer that provides insight into the mechanisms of disease initiation and progression [36]. DNA methylation is the most studied deregulated epigenetic mechanism in various types of cancer [37], especially in oral cancer [21,38]. DNA methylation regulates gene expression through different molecular mechanisms, one of which is its ability to prevent the transcription factors binding and directly influence chromatin folding [37]. To our knowledge, this represents the most extensive methylomic study in oral squamous cell carcinomas develop in PVL patients.

In our study, we did not find significant differences in either age or toxic habits between both groups of patients. Thus, we avoided these confounding factors in our results. In fact, according to the available evidence, one of the particularities of PVL is the absence of an aetiopathogenic association with tobacco [1,2,7]. Therefore, alterations in the methylation status would not be related to tobacco exposure.

In a first approximation, the unsupervised principal component analysis and the hierarchical clustering of the heatmap did not reveal clear differences in the methylome profile between the two groups of patients. This could be because both groups belonged to oral squamous cell carcinomas patients and therefore their epigenome is not as differential. In fact, only global differences would be expected when comparing cancer samples and intra-individual healthy tissue [38], cancer and inter-individual healthy controls [39], or in diseases with mutations in DNA methyltransferase enzymes that would produce a global change in the methylation profile [40]. For example, in the study by Jithesh et al. [38], they compared the methylome of 43 OSCC tissues with their adjacent paired healthy tissue and identified two clearly separated groups. Similarly, in the study by Milutin et al. [39] conducted on cytological samples in which the methylome clearly differentiated the three groups evaluated (oropharyngeal and oral carcinomas, oral lichenoid disease and healthy controls). Or the study by Simo-Riudalbas et al. [40] in patients with the rare autosomal recessive disease called immunodeficiency, centromeric instability and facial anomalies syndrome who had mutations in the DNMT3B gene and showed a global reduction in DNA methylation across all chromosomes. Furthermore, it is noteworthy that in our PCA analysis, we also did not observe defined clusters in each type of tumour samples, which may further reflect the biological heterogeneity of oral squamous cell carcinomas.

The two parallel supervised bioinformatic analyses was able to discriminate PVL-OSCC from cOSCC. From these, we identified 21 differentially methylated CpGs corresponding to 14 genes that would allow discrimination between groups. Of them, three CpGs had not been assigned to any know gene (namely cg26134913, cg11369761 and cg24104268) and the rest had been associated with genes unrelated to oral cancer. In our literature search, we found consistent previous reports regarding the genetic or epigenetic dysregulation of 12 genes (86%) in the prognosis of various types of cancer, such as cervical cancer (CASP7) [41] (BRD9) [42] (FGD5) [43], hepatocellular carcinomas (NR2E1) [44], cutaneous melanoma (ARL15) [45], lung cancer (NR2E1) [46] (AGL) [47] (ZNF777) [48] (ZNF429) [49], prostate cancer (CYP11A1) [50], bladder cancer (C18orf18) [51], renal cancer (ZNF433) [52] and leiomyosarcomas (DAXX) [53]. Similarly, the involvement of some of these genes in important cellular processes such as apoptosis (CASP7) [41], transcriptional regulation (ZNF77, ZNF433, DAXX) [54], mitochondrial degradation (PITRM1) [55] or steroid biosynthesis (CYP11A1) [50] has been described. Therefore, this dysregulation could promote cancer-association cell proliferation.

Regarding the aberrant methylation status of these genes, it has been described that BRD9 promoter methylation is associated with increased overall and progression-free survival in cervical cancer [42]. On the other hand, NR2E1 gene is silenced by hypermethylation mechanisms in hepatocellular carcinomas [44] and lung carcinomas [46]. Similarly, increased methylation of CYP11A1 has been linked to the development of prostate cancer recurrences [50]. The hypermethylation of C18orf18, also known as LINC00526, has been associated with poor prognosis and survival in bladder cancer [51]. Likewise, the hypermethylation of ZNF433 have been related to the progression of clear cell renal cell carcinoma [52] and the differential methylation and overexpression of DAXX in the progression of leiomyosarcoma [53]. FGD5 is a proangiogenic gene associated with tumour progression, whose aberrant methylation has been detected in cervical cancer carcinogenesis [43]. As can be observed, some of these methylation changes have been correlated with gene expression changes. This dysregulation is consistent with our results on methylation-mediated regulation of tumour gene expression in PVL-OSCC [24].

The differential methylation pattern of PVL-OSCC included the hypomethylation of 18 CpGs corresponding to 9 genes. Two differentially methylated probes (cg15977137 and cg24794734) were similarly epigenetically altered in both bioinformatics supervised analysis. AGL, also known as GDE, participates in glycogen degradation (glycogenolysis). Mutations in AGL have been associated with glycogen storage disease type III in different ethnicities. Its role as a tumour suppressor gene in bladder and lung cancer has also been described [56,57]. Loss of expression has been correlated with poor prognosis, aggressive tumour growth, and high mortality rates in bladder cancer, although its mechanism of action is unknown [57]. It is believed that the overexpression of the enzymes SHMT2 and HAS2 with AGL loss drives increased synthesis of glycine and hyaluronic acid, leading to enhanced glucose absorption and metabolism that promotes rapid tumour proliferation and growth [57]. On the other hand, the WRB gene, which is maternally imprinted, has been shown to exhibit a conserved imprinting pattern with methylated maternal alleles in cases of Down syndrome [58]. Additionally, a significant association between strabismus and the WRB gene polymorphism (rs2244352) has been confirmed in several cohorts [58]. Finally, it has been described that ARL15 overexpression inversely influences circulating adiponectin levels, which has been linked to a higher risk and susceptibility to type 2 diabetes, coronary heart disease, and rheumatoid arthritis [59]. Recently, a bioinformatic study has provide the first evidence that the low expression of ARL15 is associated with a favorable prognosis in cutaneous melanoma [60].

Interestingly, none of the gene promoters commonly hypermethylated in cOSCC, such as p16INK4A, p14ARF, CDKN2A, MGMT or DAPK, were identified as differentially methylated between our two groups [17,21]. This could be because the methylation levels at these promoters are similar in PVL-OSCC and cOSCC.

Based on these findings, to our knowledge, we are the first to describe the aberrant methylation in these 14 genes as prognosis biomarkers in oral squamous cell carcinomas. This could be explained by the fact that no studies have evaluated this profile in malignant tumour lesions in PVL, the heterogeneity of oral carcinomas, and the use of different platforms in assessing the methylome. The identification of methylation biomarkers in oral carcinogenesis can be used for monitoring OPMD such as PVL, to drive early detection screening strategies for oral cancer, and even enable the development of new epigenetic therapy.

We are aware of the limitations of the present study, due to the limited sample size, especially for clinical and pathological comparisons. However, it is difficult to obtain fresh frozen samples from the first malignancy of PVL. At epigenomic level, we have stablished the molecular basis for future research that focuses on functional characterization of differentially methylated genes, and now we need to validate our results in independent cohorts and multicentre studies in larger cohorts.

5. Conclusions

In this case-control study we demonstrate the molecular differences between PVL-OSCC and cOSCC, supporting the hypothesis that molecularly defined subtype classification in oral cancer could improve therapeutic development and proper patient monitoring. The identified differentially methylated genes could help us better understand the molecular carcinogenesis in PVL-OSCC, and the potential use of AGL, WRB or ARL15 as prognosis biomarkers.