RAB25 expression is epigenetically downregulated in oral and oropharyngeal squamous cell carcinoma with lymph node metastasis

ABSTRACT

Oral and oropharyngeal squamous cell carcinoma (OOSCC) have a low survival rate, mainly due to metastasis to the regional lymph nodes. For optimal treatment of these metastases, a neck dissection is required; however, inaccurate detection methods results in under- and over-treatment. New DNA prognostic methylation biomarkers might improve lymph node metastases detection. To identify epigenetically regulated genes associated with lymph node metastases, genome-wide methylation analysis was performed on 6 OOSCC with (pN+) and 6 OOSCC without (pN0) lymph node metastases and combined with a gene expression signature predictive for pN+ status in OOSCC. Selected genes were validated using an independent OOSCC cohort by immunohistochemistry and pyrosequencing, and on data retrieved from The Cancer Genome Atlas. A two-step statistical selection of differentially methylated sequences revealed 14 genes with increased methylation status and mRNA downregulation in pN+ OOSCC. RAB25, a known tumor suppressor gene, was the highest-ranking gene in the discovery set. In the validation sets, both RAB25 mRNA (P = 0.015) and protein levels (P = 0.012) were lower in pN+ OOSCC. RAB25 mRNA levels were negatively correlated with RAB25 methylation levels (P < 0.001) but RAB25 protein expression was not. Our data revealed that promoter methylation is a mechanism resulting in downregulation of RAB25 expression in pN+ OOSCC and decreased expression is associated with lymph node metastasis. Detection of RAB25 methylation might contribute to lymph node metastasis diagnosis and serve as a potential new therapeutic target in OOSCC.

KEYWORDS:

• DNA methylation
• epigenetic regulation
• head and neck cancer
• metastasis
• MethylCap-Seq
• oral squamous cell carcinoma

Introduction

Oral and oropharyngeal squamous cell carcinoma (OOSCC) are the most common subtypes of head and neck squamous cell carcinomas (HNSCC) and are characterized by an overall 5-year survival < 50%.¹ This low survival rate is greatly impacted by the presence of lymph node (LN) metatasis.² Patients with metastases in the regional lymph nodes of the neck have a 5-year survival that is half that of those who do not present regional metastases.^3,4 Therefore, for treatment decision-making, it is important to accurately detect the presence LN metastasis. Currently, diagnosis consists only of clinical examination and imaging, which are known to have low sensitivity and low specificity for LN metastasis detection.^5-8 When LN metastases are detected, a neck dissection is required, but this surgical procedure is accompanied by neck and shoulder morbidity. As a result, under- and over-treatment of OOSCC patients occurs frequently.^7,9 Currently, appropriate clinical and tumor biomarkers that predict the presence of LN metastasis are lacking.

DNA methylation is a mechanism of epigenetic modification that impacts cellular phenotypes by regulating gene expression and is known to affect carcinogenesis by altering proliferation rates and DNA repair.^10,11 As a result, DNA methylation screening has been used as a tool to predict clinical outcome and therapy response in cancer patients.^10,12 Moreover, DNA methylation of several genes has been reported to have a predictive value for nodal metastasis in HNSCC, including TWIST1,¹³ IGF2,^14,15 CDKN2A, MGMT, MLH1, and DAPK.^16,17 However, the discovery of these tumor markers has not improved clinical LN metastasis detection rate.

Recently, we have reported on the identification of new DNA methylation markers that predict LN status by MethylCap-Seq.¹⁸ The combination of enrichment of methylated DNA fragments and next generation sequencing has been established as a true genome-wide assay compared to other DNA methylation screening techniques (see ref¹⁹ for a review). Using a quantitative ranking of genomic loci by likelihood of differential methylation between OOSCC with metastasis negative LN (pN0) and OOSCC with metastasis positive LN (pN+), we identified WISP1 as a hypomethylation marker associated with pN+ OOSCC.¹⁸ In the present study, we report on a new approach tailored toward identifying potentially epigenetically downregulated genes in the metastatic OOSCC phenotype. Epigenetically downregulated genes are more suitable for opening up new clinical options, as hypermethylation can be more easily detected in an unmethylated background. In addition, methylated regions are potentially suited as therapeutic targets, thanks to the emergence of epigenetic editing and demethylating agents.¹⁹

For this purpose, we used a set of 696 genes that were previously reported to be differentially expressed between 143 pN0 and 79 pN+ OOSCC. This gene signature has a validated negative predictive power of 89% for LN metastases.^20-22 We combined the expression levels of the genes in this predictive gene signature with DNA methylation data acquired by MethylCap-Seq analysis.¹⁸ Using this approach, we identified 14 genes that were simultaneously hypermethylated and downregulated in pN+ OOSCC. In this manuscript, we report on the identification of RAB25 as the highest-ranking gene and analyze the association between expression and methylation of RAB25 and the presence of LN metastases.

Materials and methods

Patient selection

We selected treatment-naive OOSCC patients who underwent a neck dissection for primary tumor resection resulting in free resection margins upon histopathological examination at the University Medical Center Groningen (UMCG), between 1997 and 2008. Pathological revision was performed for all original hematoxylin and eosin (HE)-slides formalin-fixed, paraffin embedded (FFPE) tissue blocks. All pN0 tumors were histologically confirmed or had pN0 status with >2 y LN metastasis-free follow-up. All patient and tumor characteristics are available in Supplemental Table 1. For the immunohistochemical study, 227 OOSCC tumors were used for 5 tissue-microarrays (TMA) in triplicate, as described previously.²³ All TMA contained 7 different normal tissues that served as control. Human papilloma virus (HPV) status was tested by p16 immunohistochemistry followed by high-risk HPV PCR, as previously reported.²⁴ Out of 197 OOSCC patients, 5 were HPV16 positive. A total of 192 HPV-negative patients (102 pN0 and 90 pN+) were included for further analysis. For the MethylCap-Seq study, 6 pN+ and 6 pN0 tumors matched for age and primary tumor site were selected from the total cohort. Leukocytes were acquired from healthy women for endogenous methylation and methylation background estimation.^25,26 This study was performed in accordance with the Code of Conduct for proper secondary use of human tissue in the Netherlands (www.federa.org), and relevant institutional and national guidelines were followed.

DNA isolation

DNA isolation was performed as previously reported.¹⁸ Briefly, 2 10-μm thick FFPE sections were deparaffinized in xylene and incubated in 300 µl 1% SDS-proteinase K at 60°C overnight. DNA extraction was performed using phenol-chloroform and ethanol precipitation. The acquired DNA pellets were then washed with 70% ethanol, dissolved in 50 µl TE-4 (10 mM Tris/HCl; 0.1 mM EDTA, pH 8.0), and stored at 4°C. To check the DNA structural integrity, genomic DNA was amplified by multiplex PCR according to the BIOMED-2 protocol.²⁷ Cases with products ≥200 bp were selected for further analyses. DNA used in MethylCap-Seq was measured by Quant-iT™ PicoGreen® dsDNA Assay Kit, according to manufacturer's protocol (Invitrogen). The DNA used for pyrosequencing was measured using the Nanodrop ND-1000 Spectrophotometer (Thermo Scientific). Only samples with an absorbance ratio 260/280 nm > 1.8 were selected for further testing. The number of tumor cells required for this study was set at 60%, as estimated by HE-staining of 3-μm thick sections.

MethylCap-Seq

MethylCap-Seq analysis was performed as reported previously.¹⁸ Briefly, genome-wide methylation was assessed using 500 ng of DNA fragmented by Covaris S2 (Covaris) and obtained from 6 pN0 OOSCC, 6 pN+ OOSCC, and 2 pools of leukocytes. Methylated DNA was enriched using the methyl binding domain protein MeCP2 (MethylCap-kit, Diagenode) and followed by paired-end next generation sequencing on the Illumina GA II Sequencer (Illumina). Subsequently, the enriched, captured, and sequenced reads were mapped to the human reference genome (NCBI build 37.3) using the BOWTIE software.²⁸ Only the reads that mapped to unique loci were included. Reads that exactly overlapped with each other were excluded, as identical reads are most likely the result of amplification of the same DNA fragment. Additionally, the mapped distance between the paired-ends could not be longer than 400 bp. Finally, all the mapped reads were compared to the “Map of the Human Methylome” build 2 [http://www.biobix.be/map-of-the-human-methylome/, BIOBIX (Lab of Bioinformatics and Computational Genomics), University of Ghent, Ghent, Belgium 2014], which consists of an in-house developed summary of all experimentally assessed genomic sites of potential differential methylation [called Methylation Cores (MC)]

To identify a candidate set of genomic regions differentially methylated between pN0 and pN+ OOSCC, all MC located 2,000 bp upstream to 500 bp downstream of the transcription start site (TSS) or in the first exon of an Ensemble (v65) gene were statistically compared using R³⁵ with R-package Bayseq.²⁹ The sequencing experiment proved to be underpowered in terms of sequencing depth and number of biological replicates, precluding any definite conclusions. Therefore, we focused on the identification of the most interesting set of putatively differentially methylated regions that could be validated in a subsequent setup. This led to the following 2-step MC selection method: In the first step, the number of methylated samples was determined for both groups (pN+ and pN0). A sample was called [unmethylated] if there were no reads and [methylated] if there were one or more reads. A Fisher exact test was performed to rank the MCs for differential number of methylated samples between both groups. Ties in P-values, due to the limited number of samples, were broken by secondary ranking on log fold-change methylation between groups (average methylation incremented by 1 in both groups). In contrast to our previous quantitative ranking based on differential methylation,¹⁸ this pre-selection is unaffected by the variability of the signal in the methylated group. In the second step, the Mann-Whitney-U test was applied to the 5,000 highest-ranked MCs from the first step. MCs with a P-value < 0.05 (a total of 1,709) were retained for further consideration. Finally, only the MCs associated with genes that have an annotated function in the UniProtKB/Swiss-Prot database³⁰ were selected for further analyses.

Positive and negative predictive value for the methylation status of all MCs was calculated. For each MC all OOSCC with a read count of ≥3 reads were considered as methylated and OOSCC with a read count <3 reads were considered as unmethylated. The positive predictive value was then calculated as follows: (true positive pN+ OOSCC) / (true positive pN+ OOSCC + false positive pN0 OOSCC). The negative predictive value was calculated as: (true negative pN0 OOSCC) / (true negative pN0 OOSCC + false negative pN+ OOSCC).

Gene selection

To identify epigenetically downregulated genes in pN+ OOSCC, a validated gene signature predictive of pN-status in OOSCC (published by Hooff et al.²²) was combined with MethylCap-Seq data (Fig. 1). This gene signature is based on a diagnostic microarray consisting of 696 genes and was validated on 222 OOSCC from 8 different medical centers in The Netherlands.^20-22 Genes that were found by MethylCap-Seq to be hypermethylated in pN+ OOSCC and found to be downregulated in pN+ OOSCC by microarray were selected for further analyses.

Figure 1. Strategy to identify epigenetically downregulated genes in pN+ OSCC. On the left: published gene signatures predictive of pN-status in OSCC were used to identified significantly downregulated genes in pN+ OSCC.^20-22 On the right: MethylCap-Seq was performed on 6 pN0 OSCC and pN+ OSCC.¹⁸ All reads of MCs in gene promoter regions were ranked according to the likelihood of differential methylation and an approximate FDR. The 5,000 MCs with the lowest FDR were further tested by Mann-Whitney-U. The MC associated with genes without annotated gene functions were excluded. In the middle: the gene signature and methylation data were compared to select epigenetically regulated genes in pN+ OSCC (n = 23). From these 23 genes, epigenetically downregulated genes in pN+ OSCC were selected. Based on the amount of mRNA downregulation, statistical differences in methylation between pN0 and pN+ OSCC, and positive and negative predictive value, RAB25 was selected as the most significantly epigenetically downregulated gene in pN+ OSCC compared to pN0 OSCC.

The Cancer Genome Atlas data analysis

The Cancer Genome Atlas (TCGA) validation was performed as reported previously.¹⁸ Clinical data for all HNSCC patients (n = 423) was downloaded from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) on April 7th 2013. Patients with a tumor located in “Floor of Mouth,” “Oral Cavity,” or “Oral Tongue,” with known pathological N-status, available methylation status, and mRNA data were selected (n = 147). Patient and tumor characteristics of the selected TCGA cases are presented in Supplemental Table 1. All pathological N-statuses were dichotomized for further analyses.

For methylation analysis, level 3 methylation Illumina Infinium HumanMethylation450 (450K) data was downloaded for the previously selected oral SCC (OSCC) patients from the TCGA data portal (https://tcga-data.nci.nih.gov/tcga/) on April 7th 2013. Additional Infinium 450K probe information was acquired from the gene expression omnibus (GEO) accession number GSE42409, including, distance to TSS, associated CpG island, and chromosomal localization.⁹ All probes located up to 2,000 bp upstream and 500 bp downstream of a TSS were selected for further analyses. R (version 3.0.3),³⁵ Rstudio (RStudio, Inc.), and the Lumi package³¹ were used to convert the 450K probe β-values to M-values using the beta2m function. Subsequently, all M-values were quantile-normalized by the normalizeBetweenArrays function of R package Limma.³² Using the eBayes function of the Lumi R package, all 450K probes located 2,000 bp upstream to 500 bp downstream of the RAB25 TSS (a total of 3) were statistically compared between pN0 OSCC (n = 61) and pN+ OSCC (n = 86).³¹

For expression analysis, all mRNA expression z-scores (RNA Seq V2 RSEM) from the HNSCC TCGA “provisional cancer study” were downloaded from the cBioportal public portal (http://www.cbioportal.org/public-portal/)^33,34 on April 30th 2014 and statistically compared between pN0 and pN+ OSCC by Mann-Whitney test using R. The optimal cut-off value for RAB25 mRNA levels between pN0 and pN+ OSCC was determined to be z-score: −0.4250 by ROC-curve analysis, using SPSS version 22.0.1 (IBM). For copy number and mutation analyses, all RAB25 mutation and GISTIC data from the HNSCC TCGA “provisional cancer study” were downloaded from cBioportal public portal on April 6th 2015. TCGA survival data were incomplete and varied between testing labs and were therefore not analyzed.

Spearman rank correlations between RAB25 mRNA z-scores and normalized M-values of all RAB25 probes were calculated by the basic R function cor.test.³⁵

Putative RAB25-regulating miRNAs were identified using the miRDB database (http://mirdb.org/miRDB/) on April 6th 2015 (n = 12).³⁶ Subsequently, for all miRNAs with available data (n = 6), all RNA Seq V2 RSEM (z-score threshold of ± 2), mutation, and gene copy number data for the miRNA and RAB25 were downloaded from the cBioportal public portal (http://www.cbioportal.org/public-portal/)^33,34 on April 17th 2015. In total, 5 different types of gene copy number alterations were distinguished: -2 (homozygous deletion), -1 (hemizygous deletion), 0 (no gene copy number alteration), 1 (gain), and 2 (high-level amplification).

Bisulfite pyrosequencing

Extracted genomic DNA (1 µg/sample) was sodium bisulfite-treated using the EZ DNA methylation kit (Zymo, BaseClear, Leiden, The Netherlands) according to the manufacturer's protocol. RAB25 bisulfite pyrosequencing PCR and sequencing primers were designed using Pyromark Assay design version 2.0.1.15 (Qiagen). All primer sequences and PCR conditions are available in Supplemental Table 2. Bisulfite treated DNA was amplified using the Pyromark PCR kit according to the company protocol (Qiagen). Each reaction was performed with 12.5 µl 2x PCR master mix, 200 nmol forward primer, and 200 nmol reverse primer. PCR was performed as follows: 15 min 95°C, 50 cycles of (30 sec 94°C, 30 sec 59°C, 30 sec 72°C), and 10 min 72°C. PCR products were checked on a 2% agarose gel containing 15 µl ethidium bromide. Biotinylated PCR product (15 µl ) was captured using 1 µl Streptavidin-coated Sepharose High Performance beads (GE Healthcare). Captured amplicons were then purified using the Q24 Vacuum Workstation (Qiagen), according to the manufacturer's protocol, washed with 70% alcohol, denatured using PyroMark Denaturation Solution (Qiagen), and washed with PyroMark Wash Buffer (Qiagen). The purified PCR product was then added to 25 µl 0.3 µM RAB25 sequence primers and followed by bisulfite pyrosequencing analysis using the Pyromark Q24 (Qiagen). Pyrosequencing results were analyzed using the provided Pyromark Q24 software version 2.0.6 (Qiagen). Each pyrosequencing run included 3 control samples: leukocyte DNA from healthy patients as control for normal/endogenous methylation levels; in vitro methylated (SssI digested) leukocyte DNA as hypermethylation control; and whole-genome amplified (WGA) leukocyte DNA, amplified using the Illustra Ready-To-Go GenomiPhi HY DNA Amplification Kit (GE Healthcare), as a control for unmethylated DNA.

Immunohistochemistry

FFPE tumor tissue sections (3-μm thick) were deparaffinized in xylol and rehydrated using decreasing ethanol concentrations (100%, 96%, 80%, 70%, and 50%). Antigen retrieval was performed using a citrate buffer (10 mM citric acid, 0.05% Tween 20, pH 6.0) and heated in a microwave oven for 15 min at 300 W. Endogenous peroxidase was blocked with a 0.3% H₂O₂ solution for 30 min at room temperature, followed by incubation with a mouse monoclonal antibody to human RAB25 clone 3F12F3 (Santa Cruz), diluted 1:50 in PBS with 1% bovine serum albumin, overnight at 4°C. Subsequently, primary antibody detection was achieved by incubation with Envision+ (Dako) horseradish peroxidase for 30 min at room temperature and developed with 3,3-diaminobenzidine solution (Dako) containing 0.03% H₂O₂ and counterstained with hematoxylin for 2 min. Mammary epithelial cells were used as a control for positive RAB25 expression.³⁷ The percentage of positive tumor cells was scored as previously reported.^38,39 Three RAB25 immunoreactivity intensities were also recorded: 0 (no staining); 1 (moderate staining); and 2 (strong staining). The level of intensity of the staining was scored independently by 2 blinded observers (MJAMC and MFM). Discordant results were discussed until consensus was reached or decided by an experienced HNSCC pathologist (BvdV). The optimal cut-off between high and low RAB25 positive neoplastic cells, defined as the percentage of neoplastic cells with tumors with any level of expression (moderate/strong), was determined by ROC curve analysis relative to pN-status and established as 33% RAB25 positive tumor cells. A total of 178 out of the 192 HPV-negative HNSCC were used for RAB25 immunoreactivity analysis.

Statistical analysis

Statistical analysis was performed using SPSS (IBM) and R (version 3.0.3). Associations between RAB25 expression and clinico-pathological characteristics were tested using the χ2 test. Survival was defined as the number of days between the first treatment and disease specific death (DSS) or disease recurrence (DFS) and analyzed by Kaplan-Meier curves and log rank test. All tests were performed as 2-tailed and a P-value <0.05 was considered statistically significant.

Results

RAB25 is the highest-ranking differentially methylated and expressed gene in pN+ OOSCC

To identify genes whose expression is regulated by methylation, a validated gene expression signature and methylation data were combined using a stepwise selection approach as outlined in Fig. 1. After combining the gene signature and methylation data, 23 genes were found to be present in both the differentially methylated gene panel and the differentially expressed gene panel (Supplemental Table 3).

Out of these 23 potentially epigenetically regulated genes, 20 genes were hypermethylated in the pN+ OOSCC of the UMCG panel, as identified by MethylCap-Seq. Finally, 14 of these 20 genes (ACTA1, BRUNOL4, COBLL1, GFRA1, H2AFY, IL22RA1, KRT17, LAMP3, MALL, MAST4, NDUFA10, RAB25, S100A9, and WDR13) showed both promoter hypermethylation as well as expression downregulation in pN+ OOSCC (Table 1). Of these 14 genes, RAB25 showed the highest downregulation of expression and concomitant highest rate of hypermethylation in pN+ OOSCC (Table 1). Moreover, the RAB25 read count distribution between pN0 and pN+ OOSCC showed the highest positive and negative predictive value for pN-status (Table 1 and Supplemental Table 3). Therefore, RAB25 was studied in more detail as an epigenetically downregulated gene in pN+ OOSCC.

Table 1. Epigenetically downregulated genes in pN+ OSCC. All 14 potentially epigenetically downregulated genes in pN+ OSCC compared to pN0 OSCC after cross-reference of expression microarray and MethylCap-Seq data (see Fig. 1). The positive and negative predictive value of the reads for pN+ status, associated hypermethylation, read distribution between pN0 and pN+ OSCC, and predictive value of the methylation data are illustrated. P-value for the differential DNA methylation was calculated using the Mann-Whitney-U test. Positive and negative predictive value for the methylation status of all MCs were calculated as follows: OOSCC with a read count of ≥ 3 reads were considered true positives and OOSCC with a count read<3 were considered true negatives. Subsequently, the positive predictive value was then calculated as: (true positive pN+ OOSCC) / (true positive pN+ OOSCC + false positive pN0 OOSCC). Finally, the negative predictive value was calculated as: (true negative pN0 OOSCC) / (true negative pN0 OOSCC + false negative pN+ OOSCC).

Methylation Core data				DNA Methylation data				mRNA data	Epigenetic regulation
Gene name	Chr	Distance TSS (bp)	length (bp)	P-Value	Hypermeth in	Pos. Pred.	Neg Pred.	expression pN+	mRNA & Meth
RAB25	1	−108	233	0.02	pN+	100%	86%	−0.15	↓ Negative
COBLL1	2	−1247	191	0.02	pN+	100%	75%	−0.14	↓ Negative
GFRA1	10	−809	120	0.04	pN+	100%	67%	−0.11	↓ Negative
S100A9	1	490	125	0.04	pN+	100%	60%	−0.1	↓ Negative
LAMP3	3	0	284	0.05	pN+	80%	71%	−0.09	↓ Negative
ACTA1	1	−756	273	0.01	pN+	100%	67%	−0.08	↓ Negative
KRT17	17	−296	1	0.02	pN+	100%	55%	−0.08	↓ Negative
MAST4	5	−271	57	0.03	pN+	0%	50%	−0.06	↓ Negative
IL22RA1	1	114	229	0.05	pN+	75%	63%	−0.04	↓ Negative
BRUNOL4	18	−1543	24	0.03	pN+	100%	67%	−0.03	↓ Negative
NDUFA10	2	−1155	9	0.02	pN+	100%	55%	−0.01	↓ Negative
MALL	2	413	152	0.03	pN+	100%	55%	−0.01	↓ Negative
WDR13	X	0	54	0.05	pN+	100%	55%	−0.01	↓ Negative
H2AFY	5	−1065	90	0.03	pN+	100%	55%	−0.01	↓ Negative

Validation of epigenetic regulation of RAB25 in the independent TCGA cohort

Our data revealed a strong association between decreased mRNA expression and increased methylation of the RAB25 gene in pN+ OOSCC compared to pN0 OOSCC. To confirm this association, we selected all 147 OSCC available in the public TCGA database with available RAB25 mRNA levels, RAB25 methylation data, and pN-status data. Among the Illumina Infinium 450K probes, 5 probes were associated with the RAB25 gene (Supplemental Table 4). In total, 3 probes (cg15896939, cg09243900, and cg19580810) were located in the RAB25 promoter region (Supplemental Fig. 1). Methylation status of these 3 RAB25 promoter probes (cg15896939, P = 0.003; cg09243900, P = 0.023; and cg19580810, P < 0.001) was significantly higher in the OSCC with low RAB25 mRNA levels (Fig. 2A). Additionally, methylation levels of all 3 RAB25 probes showed a significant negative correlation with RAB25 mRNA levels (cg15896939, R = −0.230, P = 0.005; cg09243900, R = −0.162, P = 0.049; and cg19580810, R = −0.390, P < 0.001; Fig. 2B). Analysis of TGCA database confirmed that methylation of RAB25 is associated with decreased expression levels. Additionally, the location of 2 of these 3 probes (cg15896939 and cg09243900) overlapped with the RAB25 MCs annotated by MethylCap-Seq (Supplemental Fig. 1).

Figure 2. RAB25 mRNA levels in relation with the 3 RAB25 TSS 450K probes (cg09243900, cg15896939, and cg19580810) methylation levels in the TCGA OSCC cohort. (A) RAB25 methylation levels compared between OSCC with high RAB25 mRNA levels and OSCC with low RAB25 mRNA levels. The M-values of the 3 RAB25 Infinium 450K promoter probes were significantly higher in OSCC with low RAB25 mRNA z-scores compared to OSCC with high RAB25 mRNA z-scores. (B) Spearman correlations between RAB25 methylation and RAB25 mRNA levels. All 3 RAB25 promoter probes showed a significant negative correlation between RAB25 promoter probe M-values and RAB25 mRNA z-scores.

Association between RAB25 methylation and lymph node status

To determine whether RAB25 promoter methylation is associated with pN-status in OSCC, we analyzed the methylation levels of the 3 RAB25 promoter probes (cg09243900, cg15896939, and cg19580810) in 61 pN0 and 86 pN+ OSCC in the TCGA database. No significantly different methylation was found for any of the 3 RAB25 promoter probes between pN0 and pN+ OSCC (Supplemental Fig. 2A). Additionally, RAB25 methylation was measured in an independent UMCG OOSCC cohort (n = 47) using 3 different bisulfite pyrosequencing assays of the promoter region containing the annotated RAB25 MCs (bisulfite primer locations are shown in Supplemental Fig. 1). No significant differences in RAB25 methylation levels were found between pN0 and pN+ OOSCC for any of the 9 CpG sites (Supplemental Fig. 2B). These data suggest that DNA methylation of RAB25 promoter region is not directly related to LN metastasis in OOSCC.

Association between RAB25 expression and lymph node status

To determine the association between RAB25 expression and LN status in OOSCC, we analyzed RAB25 mRNA levels in OSCC using data available in the public TCGA database. Analyses of RAB25 mRNA levels in 147 OSCC revealed significantly lower (P = 0.015) RAB25 expression in pN+ (n=86) than in pN0 OSCC (n=61) (Fig. 3A). High RAB25 mRNA expression was found to be significantly associated with pN0-status (P = 0.006) (Table 2A). High RAB25 mRNA expression was also associated with decreased lympho-vascular invasion (P = 0.029) (Table 2A).

Figure 3. RAB25 expression levels between pN0 and pN+ OSCC in the UMCG and TCGA OSCC cohort. (A) pN+ OSCC in the TCGA cohort (n = 86) express significantly less RAB25 mRNA than pN0 OSCC (n = 61), as revealed by Mann-Whitney-U test. (B) pN + OSCC in the UMCG cohort (n = 87) have significantly less RAB25-positive tumor cells than pN0 OSCC (n = 91), as revealed by Mann-Whitney-U test.

Table 2. Correlations between RAB25 expression and tumor characteristics. A) Associations between RAB25 mRNA expression and the clinical characteristics of the TCGA OSCC cohort. B) Associations between RAB25 protein expression and the clinical characteristics of the UMCG OSCC cohort.

	A) RAB25 mRNA in TCGA cohort			B) RAB25 protein in UMCG cohort
	Low RAB25 mRNA	High RAB25 mRNA	P-value	Low RAB25 protein	High RAB25 protein	P-value
	N (%)	N (%)		N (%)	N (%)
Total tumors	58 (40)	89 (60)		18 (10)	160 (90)
Total patients	58 (40)	89 (60)		18 (10)	160 (90)
Gender
Male	17 (36)	30 (64)	0.576	15 (14)	96 (86)	0.053
Female	41 (41)	59 (59)		3 (5)	64 (95)
Age at diagnosis (yrs)
Median	61	60	0.412	59	64	0.197
Range	26–85	19–87		38–80	25–94
Site
OSCC	[Not Applicable]			14 (9)	142 (91)	0.18
Other				4 (18)	18 (82)
pT status
1–2	24 (41)	35 (59)	0.804	12 (10)	106 (90)	0.972
3–4	34 (39)	54 (61)		6 (10)	54 (90)
pN status
0	16 (26)	45 (74)	0.006	3 (3)	88 (97)	0.002
+	42 (49)	44 (51)		15 (17)	72 (83)
Extranodal spread (only pN+)
No	19 (48)	21 (52)	0.154	9 (19)	38 (81)	0.61
Yes	17 (65)	9 (35)		6 (15)	34 (85)
Perineural invasion
No	18 (35)	33 (65)	0.289	10 (9)	106 (91)	0.595
Yes	31 (45)	38 (55)		5 (11)	39 (89)
Lymphovascular invasion
No	27 (32)	56 (68)	0.029	12 (10)	112 (90)	0.573
Yes	17 (55)	14 (45)		3 (14)	19 (86)
Histological differentiation
Well	4 (22)	14 (78)	0.11	2 (5)	38 (95)	0.181
Moderate or Poor	54 (42)	75 (58)		16 (13)	112 (87)
Infiltration depth (mm)						0.537
Median	[Not Available]			9	15
Range				3.1 – 22	0.07 – 40
Infiltration depth (mm)
<4 mm	[Not Available]			3 (11)	24 (89)	0.823
>4 mm				13 (10)	121 (90)

To validate whether RAB25 protein expression was also associated with lymph node status in our UMCG OSCC cohort, immunohistochemistry was performed on 192 HPV-negative OOSCC. We could score RAB25 immunoreactivity in 178 OOSCC. RAB25 immunohistochemistry (example in Fig. 4) revealed a significant lower number of neoplastic cells showing RAB25 protein expression in pN+ OOSCC (P = 0.012; Fig. 3B). Using a cut-off of 33% RAB25-positive neoplastic cells to define [low] and [high] expression, low RAB25 expression was significantly associated with pN+ OOSCC (P = 0.002; Table 2B). The association between low RAB25 expression and pN+ status is in good agreement with the TCGA analysis (Table 2 and Fig. 3).

Figure 4. Representative examples of RAB25 expression in 2 OSCC as detected by immunohistochemistry. Tissues were scored by the amount of RAB25-positive cells. (A) Example of a well-differentiated OSCC with a high amount of RAB25-expressing cells. (B) Example of a poorly differentiated OSCC with a very low amount of RAB25-positive cells.

RAB25 protein level of expression was not associated with other clinical characteristics (Table 2B), DSS (P = 0.232), or DFS-survival (P = 0.260). These data support an anti-invasive function of RAB25 expression in OOSCC. Analysis of RAB25 protein levels and RAB25 MC levels revealed no associations between RAB25 methylation and RAB25 protein expression in the UMCG cohort (data not shown).

Analysis of RAB25 gene copy number alterations, mutations, and miRNAs

RAB25 mRNA expression is significantly associated with the methylation status of the RAB25 promoter (Table 1 and Fig. 2). Additionally, both RAB25 mRNA and RAB25 protein expression are associated with pN-status (Fig. 3, Table 2). However, RAB25 methylation status (Supplemental Fig. 2) is not associated with pN+ status. Therefore, DNA methylation only partly explains the regulation of RAB25 protein expression. To assess the frequency of other (epi)genetic changes that might regulate RAB25 protein expression and a possible association with LN status, the frequency of RAB25 gene mutations and gene copy number alterations were assessed in 147 OSCC cases selected from the TCGA database. We found a single OSCC case (1/147) carrying a RAB25 mutation (RAB25-Q98H). We observed RAB25 copy number gain in 29 OSCC cases (1 case with high level amplification) and hemizygous RAB25 deletions (none homozygous) were detected in 15 OSCC. RAB25 mRNA levels were significantly higher in OSCC with RAB25 gene copy number increase (P = 0.024), but RAB25 mRNA levels were not associated with hemizygous deletions of RAB25 (P = 0.330). Additionally, pN-status was not associated with RAB25 copy number gain (P = 0.540), RAB25 copy number loss (P = 0.785), or with RAB25 mRNA levels, RAB25 copy number gain (P = 0.143), or RAB25 copy number loss (P = 0.584).

The miRDB database contains 12 miRNAs putatively targeting RAB25 mRNA (hsa-miR-504-5p, hsa-miR-4725-5p, hsa-miR-608, hsa-miR-4651, hsa-miR-185-3p, hsa-miR-4520-3p, hsa-miR-4447, hsa-miR-8071, hsa-miR-4761-3p, hsa-miR-1296-3p, hsa-miR-6862-5p, hsa-miR-4253). For 6 of those miRNAs, expression, mutation, and copy number data were available from the TCGA database. None of the 6 miRNAs displayed aberrant gene expression, mutations, or copy number alterations in the 530 HNSCC present in the TCGA database (data not shown).

Discussion

We used a combination of genome-wide methylation analysis and a validated gene signature predictive for pN+ status in OOSCC to identify potential epigenetically regulated genes in the OOSCC metastatic phenotype. Of all analyzed genes, RAB25 is the most likely epigenetically regulated and predictive for pN+ OOSCC. RAB25 is reported to be a tumor suppressor gene lost in HNSCC subtypes,^38,39 that is also hypermethylated in HNSCC cell lines compared to healthy tissue,^38,39 underlining its importance and epigenetic inhibition of RAB25 protein expression during carcinogenesis.

The RAB25 protein is a member of the RAB11 subfamily of small GTPases. These GTPases are emerging as novel and important regulators of cancer development and progression (for a review see ref. 38). Aberrant expression of small GTPases in general, and RAB25 specifically,^40,41 has been detected in various cancers,^42,43 including HNSCC and OSCC.^38,39 Interestingly, changes in RAB25 expression are correlated with tumor invasiveness in almost all cancer types,^44-47 but only in triple-negative breast and HNSCC, RAB25 functions as a tumor suppressor gene, and loss of RAB25 leads to increased migration and invasion.^38,47-49

Epigenetic downregulation of RAB25 was reported in ovarian cancer compared to normal ovarian tissue,⁵⁰ esophageal cancer cell lines, compared to paired normal esophageal tissue,³⁸ and in HNSCC cell lines.^38,39 This supports the hypothesis that loss of RAB25 expression in pN+ OSCC is caused by hypermethylation, since both increased hypermethylation⁵¹ and metastasis are associated with progressive cancer⁵² and, specifically, HNSCC.⁵³ Additionally, epigenetic regulation of the expression of other small GTPases, to which RAB25 belong, has been shown in metastatic lung cancer⁵⁴ and in colon cancer.⁵⁵

We confirmed that loss of RAB25 protein expression correlated with the presence of LN metastasis in HNSCC and OOSCC, specifically,^38,39,48 and could be used to predict LN metastasis in OOSCC.^20-22 Because of our selection method, genes selected were not only associated with the presence of nodal metastases, but, more specifically, were downregulated in pN+ cases. These selection criteria are reflected in the fact that RAB25 mRNA and protein expression showed the best predictive value for nodal metastases in patients where expression is low (a high positive predictive value). Consequently, the negative predictive value is lower. This effect is independent of T-stage (data not shown), but should be taken into account in future studies into the clinical application of RAB25 expression. Additionally, MethylCap-Seq identified RAB25 as differentially methylated between pN0 OOSCC and pN+ OOSCC. These data suggest that RAB25 is epigenetically regulated and lost during cancer progression as a result of hypermethylation. However, we could not confirm differential methylation on a larger independent cohort using bisulfite pyrosequencing and Illumina Infinium 450K TCGA data, although we did find significant correlation between RAB25 mRNA levels and RAB25 DNA methylation levels. These data suggest that RAB25 is regulated by DNA methylation, but also potentially subjected to other forms of epigenetic regulation, such as histone modification or miRNAs. However, previous reports show no relation between histone modifications and RAB25 expression in esophageal cancer³⁸ and alterations of 6 miRNAs that regulate RAB25 were found to be almost non-existent in the TCGA OSCC database (this paper). Most RAB25 gene copy number alterations were amplifications and can thus not be responsible for downregulated RAB25 protein expression. The frequency of RAB25 loss was, however, very low in the TCGA OSCC database to draw firm conclusions.

In summary, our data suggest that epigenetic silencing of RAB25 contributes to LN metastasis in OOSCC patients. Therefore, RAB25 protein expression assessment might contribute to better patient diagnosis and RAB25 epigenetic editing might open new therapeutic options for the treatment of LN metastasis through epigenetic editing or using demethylating agents in order to increase OOSCC patient prognosis and care. Genome-wide methylation analysis using the MethylCap-Seq is a promising approach to identify important epigenetically regulated genes in carcinogenesis.

Disclosure of potential conflicts of interest

The CTMM Air Force Consortium is a private/public consortium with involvement of academia, private companies, and the government. It is not a commercial source of funding. M Clausen and M Mastik are funded by CTMM consortium to perform the research project (global methylation screening in HNSCC) that is a part of Clausen's thesis. These funds had no role in study design, data collection and analysis, decision to publish and preparation of the manuscript. No other external funding sources for this study. W van Criekinge is an employee for MDXHealth. The other authors have declared no conflicts of interest.

Funding

This work was funded by the CTMM Air Force consortium (http://www.ctmm.nl). CTMM pays for the salaries of M Clausen and M Mastik (partly), S Denil is supported by an IWT grant (SB101371).