Oral and pharyngeal cancer, grouped together, is the sixth most common cancer in the world. The annual estimated incidence is approximately 600,000 per year, two thirds of these cases occurring in developing countries [1]. The prevalence of oral cancers is high especially in South and Southeast Asia, where distinct cultural practices such as betel-quid chewing and varying patterns of tobacco and alcohol use are important risk factors that predispose to cancer of the oral cavity. The mortality rates of these tumors have remained unchanged (50% within 5 years after diagnosis) and are related to tobacco smoking and alcohol intake. Oral cancer patients are usually diagnosed at an advanced stage (two thirds are III–IV), which is associated with a worse prognosis and higher radio- and chemotherapy morbidity. Moreover, quality of life is disproportionately compromised in the oral cavity patient since surgical therapy can be mutilating and often has significant effects on swallowing, speech, and physical appearance. Evidently, improved oral cancer prevention, early detection, and better diagnostic and clinical management tools are needed to identify high-risk patients, such as those with smoking and alcohol exposure, patients without adequate access to health care, and patients with high-risk lesions such as oral leukoplakia (OL), which may progress to carcinoma. The estimated prevalence of OL is approximately 0.5% worldwide. Oral squamous cell carcinoma (OSCC) is the most common neoplastic disease in the head and neck region and is commonly preceded by oral premalignant lesions (OPML) such as OL. In Western countries, the annual malignant transformation rate from OL to OSCC varies from 0.13 to 36.4% [2, 3].

In addition, OSCC patients can develop a second primary OSCC, with a frequency ranging between 17 and 30% [4]. Clinical and histological features of OPML do not provide enough information to identify premalignant lesions at high risk of malignant transformation or patients who will develop an OSCC during follow-up. OSCC or high-grade squamous intraepithelial lesions (HGSIL) are usually diagnosed by incisional biopsy. Nevertheless, the biopsy requires a minimally invasive surgical approach that can create discomfort and may be refused by the patient. Moreover, only one third of patients present with early stage disease (T1–2, N0), and for the remainder, life expectancy is very short. Different researchers, including our group, have recently proposed an attractive strategy to reduce the burden of OSCC by analyzing the methylation status of a panel of genes starting from saliva and/or brushing specimens [5–10]. In many cancers, gene silencing by promoter methylation seems to be an early event in carcinogenesis and may occur even more frequently than structural inactivation of genes by mutations and deletions [11]. Histologically normal tissue adjacent to tumors and OPML can have an aberrant methylation pattern in candidate genes, suggesting that these epigenetic modifications are early events in oral carcinogenesis [12]. This study assessed the methylation status of a set of 18 gene promoters and long interspersed element 1 (LINE1) to detect early OSCC and OPML by bisulfite next-generation sequencing (NGS), starting from minimal invasive collection specimens obtained by oral brushing.

Screening populations for the early detection of asymptomatic carcinoma or precursor lesions is an attractive strategy to reduce the burden of OSCC. A major example of this approach is the cytological smear test in women for cervical cancer screening. However, low sensitivity and specificity have precluded the widescale adoption of microscopic cytology for the detection of oral cancer and OPML [38]. Methods such as toluidine blue staining and autofluorescence imaging have been proposed to improve the clinical identification of high-risk OPML, but a recent meta-analysis confirmed there was no evidence to support the use of these adjunctive technologies as screening tools to reduce oral cancer mortality [39]. Recent findings indicated that DNA methylation analysis of a specific set of genes may serve to detect early stage oral cancer lesions [22, 40], while the location of core regions and the density of methylation are required for gene silencing [41].

This study adopted a different approach from 450k methylation arrays to interrogate a set of genes already known to be altered in HNSCC [5, 15, 19, 20, 27]. For in-depth investigation of the true methylation level, we developed a novel accurate, sensitive, and specific assay to detect early OSCC and HGSIL from oral brushing specimens using bisulfite NGS. Quantitative DNA methylation analysis of the following 13 genes: ZAP70, ITGA4, KIF1A, PARP15, EPHX3, NTM, LRRTM1, FLI1, MIR193, LINC00599, MIR296, TERT, and GP1BB, clearly discriminated OSCC and HGSIL from healthy donors or from contralateral normal mucosa from the same patients. In particular, we found in OSCC and HGSIL hypermethylation of ZAP70, ITGA4, KIF1A, PARP15, EPHX3, NTM, LRRTM1, FLI1, MIR193, LINC00599, PAX1, and MIR137HG and hypomethylation of TERT, MIR296, and GP1BB at various CpGs. No epigenetic changes were found in DNMT1, H19, TERC, or any global hypomethylation of repetitive LINE1. Best performances were found for ZAP70, GP1BB, PAX1, LRRTM1, KIF1A, and TERT showing the best sensitivity and specificity in a ROC curve analysis (see the best CpGs for each gene in Figs. 3, 4, 5, and 6). Despite its high discriminant potential, we excluded the new panel PAX1 described by Guerrero-Preston [20] because it was patented by the same group, and MIR137HG which yielded ambiguous results (see Additional files 1 and 3) and was previously found to be differentially methylated in OLP [15, 42].

The stratification among distinct groups (OSCC vs healthy normal individuals, OSCC vs normal contralateral mucosa) was based not simply on evaluation of the methylated/unmethylated allele but also on quantitative identification of the correct threshold among various CpGs of the different genes considered. Combining the methylation level of 13 genes, the discrimination power of this new assay reached an optimal level of accuracy (AUC 0.981). The only OSCC case negative for the score was the OSCC with sarcomatoid features. No aberrant methylation pattern was found in any of the genes investigated, signifying that this rare variant of oral carcinoma does not share the same epigenetic status as the other OSCC cases, as shown in Fig. 8 in cluster analysis. This result may be considered for clinical management, and further investigation into genomic alterations is required.

Linear discriminant analysis allowed us to generate a score that weighted the best CpGs from the most informative genes investigated in the study. A single gene discriminating method or a single CpG may be inadequate for our purpose. Considering only one gene at a time such as GP1BB, the AUC varied between 0.968 and 0.876 depending on the CpGs, while for TERT, it varied from 0.924 to 0.667, implying the central role of CpG location (see Table 3 for the best and worst AUC values).

By using an independent cohort of samples, we further validate our algorithm confirming a clear negative value for all normal donors enrolled, together with one oral fibroma and 12 OLP. On the contrary, values over the threshold were detected in the two OSCC, in all three PVL, and in two OLP. Positivity for all PVL may be related to its nature of a progressive, multifocal, exophytic form of leukoplakia with high rates of malignant transformation. For two OLP positives (14.3%), the World Health Organization (WHO) has considered OLP as a premalignant lesion. The range of malignant transformation varies between 0 and 10% [43], which is consistent with our findings.

Basically, the calculated thresholds of most interrogated CpGs were proximal to 0 in the following genes: FLI1, ITGA4, KIF1A, MIR193, NTM, LINC00599, PAX1, and MIR137. Overall, the aberrant methylation pattern of these OSCC-related markers was estimated to be 0.33 ± 0.3 SD, simply overcoming the basal unmethylated condition seen in normal cells from both the contralateral mucosa of the same patients and the oral mucosa of healthy donors. By contrast, EPHX3, LRRTM1, and PARP15 showed higher values both for threshold and mean values for the OSCC group (0.61 ± 0.31 SD). This implies that the basal level of methylated CpG in these genes in normal tissues is not close to 0, and a small fraction of epialleles are methylated. Moreover, methylation levels in OSCC exceeded 50%. The behavior of a set of five genes including ZAP70, GP1BB, H19, EPHX3, and MIR193 fluctuated among the various CpGs as shown in Additional file 1. Interestingly, the gap between normal and OSCC remained the same (Kruskal-Wallis P values were mostly < 0.05) but the absolute values changed conspicuously among different positions investigated. This implies that each CpG must be considered apart, and a specific threshold calculated. Recently, van Vlodrop et al. [41] emphasized the emerging evidence on the importance of the location of CpG methylation in relation to gene expression and associations with clinicopathologic characteristics in cancer. Crucial CpG islands or single CpGs for regulating gene expression not only are often situated around the transcriptional start site but can also be observed more upstream or downstream.

We used bisulfite NGS methods because of their advantages: (i) many CpG sites, so that complex methylation patterns of individual DNA molecules can be determined; (ii) the longer reads can be aligned to the reference sequence more easily and accurately, especially in repetitive regions of the genome; and (iii) the long reads are more likely to cover more genotype information like single nucleotide polymorphisms (SNPs) in the neighborhood of cytosines for correlation analysis between DNA methylation and genotype. These features identified not simply different methylation levels among groups but for instance a different pattern of methylation between OSCC and HGSIL. Usually, OSCC revealed a homogeneous pattern with most epialleles methylated in cis, whereas HGSIL disclosed an irregular pattern of CpG methylation as shown in Fig. 2 for case 31. This behavior was shared among all five HGSIL for all the genes involved and is probably due to partial epigenetic modifications which reach a complete pattern only in full-blown OSCC. This behavior is visible only using bisulfite NGS and not methylation-sensitive PCR or qPCR, or Infinium™. A previous report by Bock et al. [44] compared the most promising assays for measuring DNA methylation in large cohorts, clinical diagnostics, and biomarker development, including amplicon bisulfite sequencing, enrichment bisulfite sequencing, Infinium™, EpiTyper®, and Pyroseq™. Best performances were obtained using amplicon bisulfite sequencing with high accuracy and robustness. In addition, this approach guarantees high throughput and a variable number of targets to interrogate, depending on the assay design [45].

We did not identify LINE1 global hypomethylation in OSCC with respect to the normal donors and normal contralateral mucosa as previously described [25] using pyrosequencing technique. With respect to the pyrosequencing approach which interrogates only four CpGs, as indicated by Ogino et al. [46], we evaluated more CpGs (21), with high number of reads per single specimen (hundreds). However, probably, we select a CpG-rich region which may not be involved in changes in methylation.

From a bioinformatic point of view, bisulfite NGS analysis is complex and the pipeline requires many steps. Commonly, NGS runs produce FASTQ files as output already trimmed for multiplex identifier or MID/IonXpress/Index (a short barcode sequence used to label samples/patients when multiplexing) to recognize loaded specimens. Firstly, these FASTQ files are processed for quality control (> Q30 for Illumina) and read lengths (> 100 bp) to discard primer dimers, then FASTQ files are converted to FASTA. With the advent of cloud computing and the availability of NGS webtools such as Galaxy Project [47], these steps have minimal PC requirements and are user friendly. Using perl to recognize specific sequence motifs, several FASTA files were created for each gene for each patient and then loaded onto a visualization tool such as BiQ Analyzer [48], MethVisual [49], QUMA [33], or BISMA [34]. However, these tools do not scale up with massively parallel sequencing having been designed for Sanger sequencing. Newer tools such as Bismark [50] and BS-Seeker [51] have been utilized more efficiently and can handle larger datasets generated by NGS technologies. In addition, a web application service named bisulfite sequencing pattern analysis tool (BSPAT) uses Bismark’s read alignments and methylation calling functionalities to provide further quality control, co-occurrence pattern analysis, simple allele-specific methylation analysis, visualization, and integration with other databases and tools [31]. BSPAT gives the correct methylation ratio of each epiallele investigated to be included in subsequent statistical analysis. We therefore used this tool to analyze our data as the best option, maintaining BISMA as a supporting tool for confirmation.

NTM, LRRTM1, MIR296, H19, and DNMT1 are known to be imprinted genes and are included in the imprinting gene database (http://www.geneimprint.com). However, no normal healthy donors, contralateral normal mucosa or DNA from whole blood of a pool of healthy female donors was either hemimethylated or fully methylated in DNMT1, LRRTM1, and NTM in our series. By contrast, H19 showed common signs of imprinting with half of the epialleles methylated in normal donors and contralateral normal mucosa. We found no evidence for an altered methylation pattern in OSCC cases with a mean methylation value of 0.56 ± 0.25 SD. For MIR296, another imprinted gene, a mean value of 0.95 ± 0.05 was discovered in normal donors, while it was found to be hypomethylated in all but two cases of OSCC (mean 0.91 ± 0.07, see Additional file 2). This is in agreement with a full methylation of both epialleles, while in lesions, we found a partial demethylation of only nine out of 15 CpGs tested (Additional file 1). Aberrant expression of MIR296 was previously related to gastric [52], bladder [53], and lung cancer [54]. The same behavior was found for TERT which was fully methylated in both epialleles in normal samples and partially hypomethylated in OSCC/HGSIL. Hypomethylation of this gene was previously reported in glioblastoma [55, 56]. Five out of six CpGs of the TERT gene were informative, even if the discriminating gap among different classes was found to be very slight, but calculating the ROC curve, the AUC obtained was highly consistent (0.921, see Additional file 3 for ROC between OSCC and normal donors).

Among the other 11 genes included in our algorithm, Marsit et al. [19] and our group previously identified three (GP1BB, ZAP70, KIF1A) as possible markers [15]. The present study provided evidence for a complete demonstration of their real value with the best AUC in the ROC analysis.

GP1BB encodes heterodimeric transmembrane protein that constitutes the receptor for von Willebrand factor and mediates platelet adhesion in the arterial circulation. Mutations of GP1BB are associated with Bernard-Soulier syndrome, an extremely rare inherited bleeding disorder [57]. ZAP70 gene encodes the ζ-chain associated protein kinase 70 kDa, which is a tyrosine kinase normally expressed by natural killer cells and T cells. Hypermethylation of ZAP70 gene predicted an unfavorable disease course in terms of disease progression and overall survival in chronic lymphocytic leukemia [58]. KIF1A (kinesin family member 1A) encodes a protein that is microtubule-dependent molecular motor involved in important intracellular functions such as organelle transport and cell division [59].

The other eight markers identified here included PAX1 and PARP15 already described by Guerrero-Preston et al. [20] as related to oral cancer, FLI1 involved in Ewing sarcoma [60], rectal cancer [61], and gastric cancer [62], NTM in prostate [63], EPHX3 in salivary gland adenoid cystic carcinoma [64], ITGA4 in colorectal cancer [65] ,and MIR193 in gastric cancer [66].

A second interesting finding of the present study was the behavior of normal distant mucosa from OSCC patients. Kruskal-Wallis test and multiple range test showed that the mean methylation value obtained in the group of normal distant mucosa of OSCC patients was significantly (P < 0.01) lower than the methylation value of the OSCC area but, at the same time, significantly higher (P < 0.01) than the methylation value of normal mucosa from healthy donors. Furthermore, three out of 30 normal distant mucosa (10%) exceeded the threshold, and this behavior in our opinion may not be due to a possible contamination of few cancer cells during the collection of normal mucosa, since our algorithm is based on quantitative methylation analysis and not simply on methylated/unmethylated status.

Genetic and epigenetic changes are also detected in histopathologically clean resection fields and could cause local relapse in mucosa primarily free of cancer cells. This has been explained by Slaughter’s model of “field cancerization” [67], whereas Braakhuis et al. [68] proposed a “patch-field” model in which a stem cell acquires genetic and epigenetic alterations in the initial phase, forming a clonal unit of altered daughter cells called a “patch.” A patch will transform into an expanding field acquiring additional genetic alterations, and by virtue of its growth advantage, a proliferating field gradually displaces the normal mucosa. Assuming that mucosa is predisposed to carcinogenesis due to exposure to exogenous genotoxins, an important clinical implication is that fields often remain after surgery of the primary tumor and may lead to new cancers clinicians currently refer to as “a second primary tumor” or “local recurrence,” depending on the exact site and time interval. The areas at highest risk for development of a second squamous cell carcinoma are large, sometimes extending to the lung [69]. Genetically altered cells may escape macroscopic or histopathological examination and may require sophisticated biomolecular approaches. An altered pattern of gene methylation in morphologically normal mucosa in OSCC patients may indicate field cancerization, so that these patients could have a higher risk of developing a second primary tumor or local recurrence. Further study on the methylation pattern in surgically resected patients may expand the potential of our new assay even for prognostic applications.

In conclusion, the present study confirmed the role of epigenetic alterations and aberrant methylation DNA status in HGSIL/OSCC and also revealed an altered methylation pattern in normal mucosa distant from the OSCC area. Early diagnosis of OSCC may be important for clinical management, particularly in high-risk populations, and our novel assay based on quantitative bisulfite NGS analysis could be a highly sensitive and specific method to detect early OSCC starting from non-invasive, easy-to-perform brush sampling. Further studies with a larger cohort including more HGSIL, low-grade SIL, and OLP with 5 years of follow-up are needed to elucidate the intrinsic prognostic potential of our assay.