Immunological and classical subtypes of oral premalignant lesions.

Oral squamous cell carcinoma (OSCC) is a major cause of cancer-associated morbidity and mortality and may develop from oral premalignant lesions (OPL). An improved molecular classification of OPL may help refining prevention strategies. We identified two main OPL gene-expression subtypes, named immunological and classical, in 86 OPL (discovery dataset). A gene expression-based score was then developed to classify OPL samples from three independent datasets, including 17 (GSE30784),13 (GSE10174) and 15 (GSE85195) OPLs, into either one of the two gene-expression subtypes. Using the single sample gene set enrichment analysis, enrichment scores for immune-related pathways were different between the two OPL subtypes. In OPL from the discovery set, loss of heterozygosities (LOH) at 3p14, 17p13, TP53, 9p21 and 8p22 and miRNA gene expression profiles were analyzed. Deconvolution of the immune infiltrate was performed using the Microenvironment Cell Populations-counter tool. A multivariate analysis revealed that decreased miRNA-142-5p expression (P = 0.0484) and lower T-cell, monocytic and myeloid dendritic cells (MDC) immune infiltration (T-cells, P = 0.0196; CD8 T cells, P = 0.0129; MDC, P = 0.0481; and monocytes, P = 0.0212) were associated with oral cancer development in the immunological subtype only. In contrast, LOH at 3p14 (P = 0.0241), 17p13 (P = 0.0348) and TP53 (P = 0.004) were associated with oral cancer development in the classical subtype only. In conclusion, we identified 2 subtypes of OPLs, namely immune and classical, which may benefit from different and specific personalized prevention interventions.

Introduction

Head and neck squamous cell carcinoma (HNSCC) is ranked as the seventh most common cancer worldwide, with more than 600,000 cases diagnosed annually, and is a major cause of cancer-associated morbidity and mortality.[1] Patients with squamous cell carcinomas of the oral cavity (OSCC), the most common anatomical site of HNSCC, have the greatest risk of death,[2] which is commonly related to local-regional advanced disease.[3] Since OSCC may develop from oral premalignant lesions (OPL), prevention of malignant transformation of OPL to OSCC may substantially improve the clinical outcome of patients. Nonetheless, one of the greatest challenges in managing patients with oral OPL is determining risk for malignant transformation and the need for intervention, with a risk of OSCC between 0.13% and 34%[4] during periods of up to 30 years. The current histologic classification of oral epithelial dysplasia is subjective and unreliable, and surgical resection of OPL has not been shown to reduce OSCC development.[5]

Some biomarkers of risk of malignant transformation, such as loss of heterozygosity (LOH),[6-8] and epithelial growth factor receptor (EGFR) gene copy number gain,[9] have been identified and validated in OPL. Retrospective and prospective studies completed initially by our group[6] and subsequently by others[7,8] demonstrate that LOH profiles (primarily at 3p14 and/or 9p21) in OPL are the most robust marker of cancer risk in this setting. Despite these efforts, randomized controlled trials (RCT) of chemopreventive agents have not yet resulted in the development of an intervention that can be considered as standard of care.[10-14]

The recent Erlotinib for Prevention of Oral Cancer (EPOC)[15] trial was the first RCT to use molecular selection (LOH) to treat patients suffering from high-risk OPL. The rationale for testing erlotinib for prevention of OSCC was partly based on the significant decrease of premalignant and neoplastic lesions induced by erlotinib treatment in the 4-nitroquinoline-1-oxide (4-NQO)-induced murine model of oral carcinogenesis.[16] Although EPOC paved the way for the beginning of an era of personalized-based prevention strategies,[17] it failed to reach the primary endpoint of oral cancer-free survival.[15]

To provide insight into the molecular heterogeneity of HNSCC, four distinct gene expression-based molecular subtypes have been identified, reproduced and well-characterized in terms of their genomic alterations, canonical pathways and methylation profiles.[18-21] Interestingly, some of these molecular subtypes present different drug sensitivity patterns[21] as well as different levels of radioresistance,[22] suggesting their potential relevance to define personalized-based therapy in head neck cancers.

Using a similar rationale, we hypothesized that the identification of OPL gene expression-based molecular subtypes may help design precision prevention strategies. Using a discovery set of 86 OPL samples, we identified two distinct gene expression subtypes of OPL, namely immunological and classical, that were validated in three independent datasets. Deconvolution scores of immune populations were computed in all datasets. Loss of heterozygosity at specific chromosomal loci as well as miRNA expression were measured in the discovery dataset and their association with oral cancer free-survival was tested.

Results

Discovery of two molecular subtypes of OPL: ‘immunological’ and ‘classical’

To identify gene expression based OPL subtypes, we used genome-wide expression profiles that we previously generated in a set of 86 oral leukoplakias[23] (GSE26549). Using the 2,500 most variable genes, we performed an unsupervised and unbiased clustering of the 86 OPL from the discovery dataset, as previously described.[18] Two main clusters, including 42 (cluster 1) and 44 (cluster 2) OPLs, respectively, were identified based on the heatmap of the consensus matrix (Figure 1A) and plots (Supplementary Figure 1) produced by the ConsensusClusterPlus R package, as well as on the heatmap of hierarchical clustering using ArrayStudio software (Supplementary Figure 2).

Figure 1.

Identification of two gene expression-based subtypes, “classical” and “immunological”, in the discovery dataset of 86 OPL. To identify gene expression subtypes, we ran the ConsensusClusterPlus R package[55] with the 2,500 most variable genes, which were selected using the median absolute deviation, as previously described [18]. (A) Heatmap of the consensus matrix for k = 2 clusters, with samples in rows and columns, and consensus values ranging from 0 (never clustered together) to 1 (always clustered together) and marked by white to dark blue. (B) Enrichment scores of 1,329 canonical pathways in each OPL samples were computed and the 15 pathways most enriched in the immunological and the classical subtypes are shown. (C) The percentage of stromal inflammatory cells (SIC) as well as intraepithelial inflammatory cells (EIC) were compared between the two subtypes (Fisher’s exact test). Abbreviation: NS: not significant.

To gain some insight into the biological differences between these two clusters, the enrichment score (ES) of 1,329 canonical pathways was computed for each OPL sample, using the single sample Gene Set Enrichment Analysis (ssGSEA) tool.[24-26] The ES was then compared between the two clusters using a Wilcoxon test. Since we observed a strong and significant enrichment of immune pathways, especially lymphocyte-related pathways, in cluster 1 compared to cluster 2 (Figure 1B; supplementary file 1), cluster 1 was named ‘immunological’. Conversely, cluster 2 was characterized by a moderate enrichment of pathways involved in xenobiotic metabolism as well as EGFR signaling and was therefore named ‘classical’ by analogy with the classical subtype of head neck cancers.[18,21] Of note, EGFR gene expression was significantly higher in the classical subtype compared to the immunological subtype of OPL (P < 0.0001; Supplementary Figure 3).

As the main biological difference between the two subtypes lied in the tumor microenvironment, we used the MCPcounter R package[27] to deconvoluate the immune infiltrate and compute a score of the absolute abundance of eight immune and two stromal cell populations. When compared to the classical subtype, the immunological subtype had higher score of T-cells (P < 0.0001), cytotoxic lymphocytes (P < 0.0001), B lineage (P = 0.0025), monocytic lineage (P < 0.0001) and myeloid dendritic cells (P < 0.0001) (Supplementary Table S1). No significant differences were found when we evaluated the abundance of other cell populations such as neutrophils, endothelial cells and fibroblasts. Interestingly, the presence of intraepithelial inflammatory cells (EIC) and stromal inflammatory cells (SIC) as defined previously[28] was significantly higher in the immunological subtype compared to the classical subtype (P = 0.0014) (Figure 1C). No significant difference was observed with regard to clinical and pathological characteristics between the two subtypes (Table 1).

Table 1.

Characteristics of the 86 patients from the OPL discovery dataset. Characteristics were compared between the immunological (N = 42) and the classical (N = 44) subtypes, using a fisher’s exact test. RP: retinyl palmitate; BC: beta-carotene; 13cRA: 13-cis retinoic acid.

Variable	ALLN = 86	Immunological N = 42	ClassicalN = 44	P–value
Oral cancer dvt
Yes	35 (41)	21 (50)	14 (32)	0.1241
No	51 (59)	21 (50)	30 (68)
Gender
Male	45 (52)	20 (48)	25 (57)	0.5174
Female	41 (48)	22 (52)	19 (43)
Ethnicity
White	78 (91)	38 (90)	40 (91)	1
Other	8 (9)	4 (10)	4 (9)
Alcohol
Current	49 (57)	20 (48)	29 (66)	0.1409
Former	8 (9)	6 (14)	2 (5)
Never	29 (34)	16 (38)	13 (30)
Smoking
Current	22 (25)	11 (26)	11 (25)	0.5847
Former	35 (41)	19 (45)	16 (36)
Never	29 (34)	12 (29)	17 (39)
Age
Median	57.5	59.5	54.5	0.0949
Range	23–90	24–90	23–80
Treatment arm
BC+ RP	21 (24)	12 (29)	9 (20)	0.6716
13cRA	47 (55)	22 (52)	25 (57)
RP only	18 (21)	8 (19)	10 (23)
Histology at baseline
Hyperplasia	54 (63)	26 (62)	28 (64)	1
Dysplasia	32 (37)	16 (38)	16 (36)

Overall, we identified two biologically distinct gene expression subtypes of OPL that harbor major differences in terms of immune microenvironment, despite being clinically similar.

Classification of OPL samples using a 400-gene expression based score

To build and test a classifier of OPL samples, the discovery dataset of 86 OPL samples was randomly split into two subsets of 57 samples (training set) and 29 OPL samples (testing set) (see methods). In the training subset, we used a method of classification based on a gene-expression ES of genes which were differentially expressed between the two subtypes. The comparison of genome-wide expression profiles resulted in the identification of 200 most overexpressed genes (UP subset) and 200 most under-expressed genes (DOWN subset) in the immunological subtype compared to the classical subtype (Wilcoxon test with a Q-value< 0.05 and a |FC| > 1.5; supplementary file 2). Using ssGSEA, we computed an ES of the UP and DOWN subsets and a combined score equal to the ESUP – ESDOWN. Using this “ES-based classifier”, an OPL was classified as being immunological or classical if the score was > 0 or < 0, respectively.

In order to compare the ES-based classifier to alternative methods of classification, three well-established algorithms of classification based on centroid (partial least square, linear discriminant analysis and k-nearest neighbor) were also applied to the 57 OPL samples from the training set (see method) using 500 selected genes (F-test) and leave one-out cross validation. The best classifier in the training set was the Partial Least Square with K = 3 (i.e., PLS3 classifier) (Supplementary Figure 4). The ES-based classifier score and the PLS3 classifier were applied to the remaining 29 OPL samples from the testing set and the percentage of classification accuracy was compared between the two classifiers. The ES-based score had a higher accuracy (96.5%) compared to PLS3 (89.7%), misclassifying only 1 of the 29 OPLs from the testing subset.

Subsequently, we used the ES-based classifier in order to classify 17, 13 and 15 new OPL samples from two independent validation sets (V1,V2 and V3 respectively) which were downloaded from Gene Expression Omnibus (see method). The distribution of samples in each molecular subtype for all datasets is provided in Supplementary Fig. 5. The ES of 1,329 canonical pathways was computed and compared between the two subtypes in V1, V2 and V3 (Supplementary File 1). We computed the estimate ([Mean pathway in immunological subtype] – [mean pathway in classical subtype]) in V1,V2 and V3, and tested its correlation with the estimate calculated in the discovery dataset (Figure 2A). Coefficients of correlation were high (r > 0.6, P < 0.0001) for both V1, V2 and V3, emphasizing the consistency of the biological differences between the two molecular subtypes across multiple datasets.

Figure 2.

Validation of biological pathways and immune populations infiltrate characterizing the two subtypes. (A) In each dataset (discovery set and three validation sets: V1 = GSE30784, V2 = GSE10174 and V3 = GSE85195), the enrichment scores (ES) of 1,329 canonical pathways and the estimate was calculated for each pathway in the following manner: (Mean score Pathway in immunological score) – (Mean score Pathway in classical subtype). We tested the correlation of estimates calculated for each pathway in each validation set with the estimate of the corresponding pathways computed in the discovery set (Pearson’s correlation). (B) Using MCPcounter to deconvoluate the immune infiltrate, the score of T-cells (B) and monocytes (C) scores were computed and compared between the two subtypes in the discovery and validation sets (Wilcoxon test).

Deconvolution of the immune infiltrate was also performed in V1 and V2 and the scores of the immune and stromal populations were compared between the two subtypes using a Wilcoxon test. The T-cell and monocyte scores were higher in the immunological subtype in the two validation datasets V1 (P = 0.0553 and P = 0.0068 respectively), V2 (P = 0.414 and P = 0.0503 respectively) and V3 (P = 0.1377 and P = 0.0015 respectively) (Figure 2B–C) (Supplementary Table S1).

Overall, we defined and validated an accurate method of classification based on a gene-expression enrichment score (ES) of 400 genes that were differentially expressed between the two subtypes.

Integrated network of mRNAs and mrnas in OPL subtypes

The main difference between the two subtypes of OPL lied in the immune microenvironment. Because microRNA-mediated RNA interference has been described as an important regulatory mechanism of the immune response,[29] we analyzed miRNA expression profiles of 71/86 OPL samples from the discovery dataset (37 “classical” and 34 “immunological”). A total of 46 and 15 miRNAs were found to be overexpressed and underexpressed respectively in the immunological subtype compared to the classical subtype (Wilcoxon test, Q-value< 0.05, |FC|> 2) (Supplementary File 3). Of note, 14/46 miRNAs overexpressed in the immunological subtype, have been involved in the activation or differentiation of lymphocytes and/or monocytes.[30,31] Using the miRNet tool, we identified 3,251 validated mRNAs targets of the 31/61 differentially expressed miRNAs (Figure 3). Interestingly, 970/3,251 (30%) mRNAs targets were also found to be differentially expressed between the two classes of OPL (Wilcoxon test, Q-value< 0.05). An integrated network of these differentially expressed mRNAs and miRNAs was performed using the Gephi software (Figure 3A). Over-expressed miRNAs and their underexpressed target genes (Figure 3B), as well as the underexpressed miRNAs and their over-expressed target genes (Figure 3C) in the immunological versus classical subtypes, were connected within two distinct networks. These networks highlighted 7 miRNAs connected with over 30 target genes. Two and five miRNAs were under-expressed (hsa-mir-30c-1-3p; mir-508-5p) and over-expressed (mir-142-3p; mir-512-3p; mir-519b-3p; mir-223-5p; mir-142-5p) in the immunological versus classical subtypes respectively (Figure3 B–C).

Figure 3.

Integrative network of miRNAs and target genes differentially expressed between the two subtypes. In the discovery set, we identified a set of 31 miRNAs differentially expressed between the immunological and classical subtypes (Q-value< 0.05, |FC|> 2) and connected with 970 target genes (miRNet tool) which were also differentially expressed between the two subtypes (Q-value< 0.05). The Gephi software allowed visualizing the integrative network of these 31 miRNAs and their 970 target genes (A). Two central sub-networks were identified: (B) over-expressed miRNA connected to under-expressed targeted genes; (C) under-expressed miRNA connected to over-expressed targeted genes. In this network, miRNAs are shown as purple circles with increasing size according to the number of connections with its target gene (yellow circles). Connections between a miRNA and a target gene are shown as colored curved lines: purple lines = miRNA overexpressed – gene overexpressed; turquoise lines = miRNA underexpressed – gene underexpressed; red lines: miRNA overexpressed – gene underexpressed; green lines: miRNA underexpressed – gene overexpressed; yellow circles = genes; purple circles = miRNA. Nodes (miRNA and genes) were clustered using the ForceAtlas2 algorithm.

In summary, the two gene expression-based subtypes of OPL were also characterized by distinct patterns of miRNAs, some of which known to be involved in the regulation of the immune response.

Molecular subtypes of OPL and oral cancer risk

We then tested whether the immunological and classical subtypes of OPL were clinically relevant using the discovery cohort in which samples were collected in the context of a clinical trial.[13] No difference was observed in OCFS between patients with immunological versus classical OPL (Supplementary Fig. 6).

We then looked at effect of LOH at specific chromosomal loci, a previously reported biomarker of risk of oral cancer development in patients with OPL.[6,7,15,32] LOH at 3p14, 3p21, 17p13, TP53 and 8p22 were analyzed as previously described in 70/86 (81.4%) OPL samples from the discovery cohort.[6] The percentage of LOH at these sites was not significantly different between the two subtypes (P > 0.05, Fisher’s exact test) (Table 2). We tested the association of LOH at these sites with OCFS, using Kaplan Meier curves (Fig; 4A-B) as well as a univariate and multivariate cox proportional hazard model including age, treatment arm and histology (hyperplasia vs. dysplasia) (Supplementary File 4). No significant association was observed in the overall population of patients with OPL. However, LOH at D17S1176 (17p13) (P = 0.0348; Odds Ratio (OR) = 12.40; Confidence Interval 95% (CI95) = [1.197;128.497]), at TP53 (P = 0.004; Odds Ratio (OR) = 52.19;CI95 = [3.612;754]) and D3S1285 (3p14) (P = 0.0241; OR = 6.46; CI95 = [1.277;32.681]) were associated with worse OCFS in the classical but not in the immunological subtypes.

Table 2.

Percentage of loss of heterozygosities (LOH) in the discovery dataset. LOH were evaluated in the discovery dataset (n = 70) at different chromosomal loci: D9s171 and D9S1747 (9p21), D3s1285 (3p14), D17s1176 (17p13), TP53 and D8s254 (8p22). Percentage of LOH was compared between immunological and classical OPL (Fisher’s exact test). “Missing data” actually include both non-informative cases and cases for which analysis could not be performed.

Loci	ALLLOH	ALLMissing data	ImmunoLOH	classicalLOH	p-value(immuno vs classical)
D9s171	12 (17.1%)	34 (48.6%)	6 (17.1%)	6 (17.1%)	1
D3s1285	14 (20.0%)	26 (37.1%)	7 (20.0%)	7 (20.0%)	1
D17s1176	12 (17.1%)	29 (41.4%)	7 (20.0%)	5 (14.3%)	0.7337
D8s254	12 (17.1%)	30 (42.9%)	8 (22.9%)	4 (11.4%)	0.5048
D9s1747	9 (12.9%)	33 (47.1%)	7 (20.0%)	2 (5.7%)	0.2616
TP53	11 (15.7%)	22 (31.4%)	3 (8.6%)	8 (22.9%)	0.1733

In addition, we evaluated the association of 8 immune (T-cells; CD8 T-cells; cytotoxic lymphocytes; NK cells; B lineage; monocytic lineage; myeloid dendritic cells; and neutrophils) and 2 stromal (endothelial cells and fibroblasts) cell populations based on their deconvolution scores with OCFS, using Kaplan-Meier curves (Figure 4 C-D) and a multivariate cox proportional hazard model including age, treatment arm and histology (hyperplasia vs. dysplasia) (Supplementary File 4). While no significant association was observed between immune or stromal cell populations and OCFS in the overall population, increased scores of the monocytic lineage (P = 0.0212; OR = 0.22; CI95 = [0.06;0.80]) and myeloid dendritic cells (P = 0.0481; OR = 0.55; CI95 = [0.30;0.99]) were associated with improved OCFS in the immunological (n = 42) but not in the classical (n = 44) subtypes (Figure 4C–D and Supplementary File 4:). Similarly, increased scores of T-cells (P = 0.0196; OR = 0.3264; CI95 = [0.13;0.84]) and CD8 T-cells (P = 0.0129; OR = 0.06; CI95 = [0.007;0.552]) were associated with improved OCFS (Supplementary File 4).

Figure 4.

OPL subtypes and biomarkers of oral cancer risk in the discovery set. In the 86 OPL discovery dataset, classical (A, C, E) and immunological (B, D, F) OPL were split into two groups according to loss of heterozygosity at D17S1176 (17p13) (A-B); the median score of monocytic lineage (C-D) (high i.e. > median), the median expression of miR-142-5p (E-F) (high i.e. > median). Oral-cancer free survival curves were compared between these groups, using a log-rank test.

Finally, we tested whether the 7 miRNAs identified in the integrative analyses were associated with the OCFS. Consistent with the association of the monocytic lineage with improved OCFS in the immunological subtype only, hsa-miR-142-5p monocyte-macrophage transducer was the only one significantly associated with a lower risk to develop OSCC in the immunological subtype only (P = 0.0484; OR = 0.68; CI95 = [0.46;0.997]), while no association was found in the classical subtype (Kaplan Meier curves: Figure 4E–F; univariate and multivariate cox models: Supplementary File 4).

Discussion

Large-scale programs have provided important biological insights into HNSCC and allowed for the identification of four gene expression-based subtypes[18,19] that have been associated with different patterns of drug sensitivity[21] and response to radiation therapy.[22] A pre-cancer atlas is the next step for elucidating the molecular heterogeneity at the early steps of head and neck tumorigenesis and help designing next generation prevention interventions.[33] Given that 0.13% to 34% of OPLs will transform into invasive SCC[4] over a period of up to 30 years, the generation of prospective collections of OPL samples has been challenging. Using previously established independent cohorts of OPL, we defined a gene expression-based classification of OPL into two subtypes, immunological and classical. The immunological subtype was characterized by overexpression of miR-142-5p as well as increased infiltration of monocytes and T-cells according to computed deconvolution of immune populations, all being associated with OCFS in this subtype only but not in the classical subtype.

Molecular classifications have been reported in different cancer types, using well-established classification algorithms, such as linear discriminant analysis, partial least square and k-nearest neighbor.[34,35] In our study, we applied these algorithms to build a classifier of OPL that we compared to an alternative classification method based on a 400 gene expression score using the ssGSEA tool. Scores based on gene expression and computed in individual samples have previously been proposed for classification purposes in biological samples e.g. epithelial – vs. mesenchymal-like samples using a pancancer signature of epithelial-to-mesenchymal transition.[36] Unlike algorithm-based methods, they allow classifying new samples independently from previous ones, which is attractive in the clinical setting, especially in the context of the forthcoming clinical trials evaluating new prevention strategies. Our OPL classifier was built from a training subset of the discovery dataset, including a relative small number of samples. Moreover, the discovery dataset includes samples from patients who had been diagnosed with oral leukoplakia and who had been randomly assigned to intervention with different chemopreventive agents, that could also limit interpretation of our results. However, the association of biomarkers with oral cancer free survival was tested using multivariate analysis including the treatment arm. Moreover, biological differences between the two gene expression based subtypes were validated in three independent datasets, which greatly enhanced confidence in our results. Nevertheless, our ES-based classifier needs to be validated in larger cohorts.

We report two subtypes of OPL while 4 subtypes of HNSCC have been proposed, i.e. classical, basal, atypical and mesenchymal-like.[18,19] More recently, the mesenchymal-like subtype has been questioned as it may be reflecting the contribution of non-malignant mesenchymal cell types of the tumor micro-environment.[37] Finally, a recent pancancer study from TCGA has reported 6 immune subtypes that encompass multiple cancer types.[38] The vast majority of HNSCC fits into the wound healing (more frequent in the classical subtype) or the IFN-gamma dominant subtypes that do not completely overlap with previously described subtypes. Whether the immunological subtype of OPL precedes transformation to the IFN-gamma subtype and the classical subtype of OPL precedes transformation to the wound healing and classical subtypes is a plausible hypothesis that remains to be demonstrated through the analysis of OPL and subsequent paired oral squamous cell carcinoma. Whether driver and passenger mutation patterns are different in these different subtypes of OPL remain to be studied. Moreover, we have recently shown that the main biological and actionable difference between oral SCC from never-smoker never-drinker and smoker drinker patients lies into the immune microenvironment. It was tempting to hypothesize that this difference would also be observed in premalignancy. However, we found no significant association between the smoking and drinking habits and immunological subtyping. The small size of our cohort may explain the lack of statistical significance of this association and larger cohorts are needed to confirm this result.

In our analyses, the immunological subtype was characterized by a strong enrichment in immune-related pathways, especially pathways related to lymphocytes and monocytes functions and a higher infiltration of monocytes and T-cells as compared to the classical subtype. This result was consistent with the overexpression of 14 miRNAs in the immunological subtype that have been involved in the activation or differentiation of lymphocytes and/or monocytes.[30,31] In particular, miR-142-5p expression was higher in the immunological subtype compared to the classical subtype while its decreased expression was associated with a higher risk of oral cancer development in the immunological subtype only. A study of microRNA profiles during progression of oral leukoplakia to cancer has shown the complex deregulation of miR-142-5p expression during oral carcinogenesis. It was found to be upregulated in cancer versus normal mucosa, while its expression changes in hyperplasia, dysplasia and cancer was not consistent across all patients, suggesting that its dysregulation could be different depending on the molecular subtype of OPL.[39] In other contexts, miR142-5p has been reported to play a role into macrophages differentiation[40] and to be down-regulated in glioblastoma-infiltrating macrophages.[41] Overall, the interaction between the immune system and neoplasia reflects a fundamental principle, applicable to all organ/cell types.[42] The recent success of the use of PD-1 immune checkpoint blockade in advanced HNSCC[43-45] confirms the clinical importance of immune modulation in this context. Overall, tumors orchestrate a complex signaling network to induce immune tolerance.[46] These events can occur as early as in the OPL, as demonstrated in preclinical[47] and early clinical studies[48] that illustrate the dynamic and complex of the premalignant immune microenvironment.

The classical subtype of OPL was characterized by an overexpression of Epidermal Growth Factor Receptor (EGFR) as well as a significant enrichment by EGFR signaling pathway. Inhibition of the EGFR pathway has been proposed for oral cancer prevention, based on the significant decrease of premalignant and neoplastic lesions induced by EGFR inhibitors in the 4-nitroquinoline-1-oxide (4-NQO)-induced murine model of oral carcinogenesis. However, the recent Erlotinib for Prevention of Oral Cancer (EPOC) trial failed to reach its primary endpoint (oral cancer-free survival). In light of our results, it is tempting to hypothesize that the benefit of EGFR targeting agents could be higher in patients suffering from classical lesions, and therefore, that the proposed molecular classification may potentially help for refining oral cancer prevention strategies.

In our study, LOH status was not significantly different between the classical or the immunological OPL subtypes. Interestingly, although LOH results are limited by the number of primers used for each region, LOH at 3p14, TP53 and 17p13 were significantly associated with a higher risk of oral cancer development in the classical subtype but not in the immunological subtype. In contrast to LOH, down-regulation of miRNA-142-5p and decreased infiltration of monocytes, myeloid dendritic cells and T-cells according to computed deconvolution of immune populations, were associated with poor OCFS in the immunological subtype but not in the classical subtype. Our results suggest that beyond LOH, our gene expression-based classification of OPL may provide an opportunity to refine the molecular assessment of oral cancer risk. Consistent with previous observations in HNSCC,[49-51] we found that premalignant-infiltrating T-cells and monocytes are associated with a lower oral cancer risk, suggesting an important role of these populations during oral carcinogenesis. Scattered evidence on the biology of OPL suggests the presence of an immune stimulatory environment characterized by IFN-gamma, and other pro-inflammatory cytokines and immune cells, including NK cells, macrophages and activated CD8 + T-cells. The 4-NQO mouse model of oral carcinogenesis, which mimics the effects of tobacco and induces molecular, genetic and tumor antigen changes similar to those of human OPL and oral SCC,[52] provides an opportunity to study the mechanisms governing the evolution of pre-malignant clones and early immune interception of this oncogenic process. Recent reports on the 4-NQO model found that OPL cells themselves elicit increased cytokine production (e.g., IFN-gamma), immune cell (e.g., activated T-cell) activation,[53] and therefore serve as a basis to define a mechanism by which OPLs modulate the immune response. Importantly, these studies suggest that immune intervention may be effective in patients with OPL as they are in a significant proportion of patients with oral cancer. Recent reports showing that metformin exerts immune-mediated anti-tumor effects,[54] in addition to preventing tumorigenesis and malignant conversion of OPL to oral cancer in the 4-NQO model, we have launched an immune interception trial assessing the activity of Metformin For Oral Cancer Prevention (M4OC-Prevent) (www.clinicaltrials.gov, NCT02581137).

In conclusion, we report a gene expression-based classification of OPL, into immunological and classical subtypes, that may be applicable in the clinical setting for refining prevention strategies. Large prospective cohorts of OPL are needed to validate our results. Further characterization of the immune infiltrate is required, which may pave the road to immune intervention strategies.

Methods

Four independent sets of OPL samples

Four independent sets of OPL samples with available gene expression profiles were established, using data generated in OPL biopsies collected prospectively in a chemoprevention trial run at The University of Texas MD Anderson Cancer Center (MDACC) (discovery set), as well as data downloaded from the public data repository Gene Expression Omnibus (GEO) (validation sets V1, V2 and V3). Written informed consent was obtained for all patients included in the trials at MDACC.

The discovery set included 86 OPL biopsies collected in patients followed prospectively with a median follow-up of 7 years (95% CI [5.6–8.6]) in a chemoprevention trial with available gene expression profiles previously published in GEO (GSE26549).[13,23] Patients had been diagnosed with OPL and randomly assigned to intervention with 13-cisretinoic acid (13cRA) versus retinyl palmitate (RP) with or without b-carotene (BC). Of note, loss of heterozygosity status at 9p21 (D9S171, D9S1747), 3p14 (D3S1285), 17p13 (D17S1176), TP53, 8p22 (D8S254) was available for 70/86 OPLs. Moreover, miRNA gene expression profiling was generated in 71/86 OPLs (see below the subsection “miRNA profiling of OPL”). Finally, based on the examination of OPL biopsies under a multi-headed microscope by three observers, each OPL was defined by the presence or absence of inflammatory cells in the stroma (SIC) and intra-epithelial inflammatory cells (EIC).[28]

In order to validate our gene expression based classification, three independent gene expression validation datasets V1, V2 and V3 were downloaded from Gene Expression omnibus and included 17 oral dysplastic lesions (GSE30784),13 OPL with hyperplasia or dysplasia (GSE10174) and 15 OPL (GSE85195) respectively.

A detailed description of the different datasets is provided in Supplementary Table S2.

MicroRNA profiling of OPL in the discovery set

All steps leading to the generation of raw data using the samples from the discovery set were processed at the University of Texas M.D. Anderson Cancer Center (MDACC) Genomics Core Facility. We produced miRNA profiles in 71/86 OPL (discovery set) using available aliquots previously prepared for mRNA gene expression profiling.[23] Each sample was hybridized on TaqMan® low density arrays Cards A and B, using the manufacturer’s instructions (Applied Biosystems, Foster City, CA). All TLDA cards were run on the Applied Biosystems 7900HT real-time PCR platform (Applied Biosystems, Foster City, CA).

Loss of heterozygosity analysis

In 70/86 OPL, loss of heterozygosities (LOH) at 9p21 (D9s171, D9S1747), 3p14 (D3s1285), 17p13 (D17s1176), TP53 and 8p22 (D8s254) were evaluated as previously reported.[6]

Unsupervised discovery of OPL gene expression subtypes

Using gene expression data of the 86 OPL samples from the discovery set (GSE26549),[23] identification of the distinct molecular subtypes of OPL was performed as previously described.[18,34] Gene expression values were gene median centered and the 2,500 most variable genes, using the median absolute deviation, were selected for clustering. We then used the ConsensusClusterPlus R package to cluster the OPL samples in an unbiased and unsupervised manner.[55] This procedure was performed with 1,000 randomly selected sets of microarray samples using a sampling proportion of 80% and a distance metric equal to one minus the Pearson correlation coefficient.

Classification methods of OPL subtypes

The discovery dataset of 86 OPL samples was randomly split into a training and a testing subsets of 57 and 29 OPL samples respectively using the “sample” function in the R language. Using the ArrayStudio© software (OmicSoft©, Cary, NC), different classification algorithms (partial least square, linear discriminant analysis, and K-nearest neighbor) were applied to the training subset to build a classifier, including the 500 most significant genes based on F-test. According to an external cross validation, the best gene classifier was selected and applied to the testing subset. In addition, based on differential gene expression between OPL subtypes, we computed a gene expression score for the OPL classification, using the single sample gene set enrichment analysis (ssGSEA)[24,25] from Gene Pattern.[26] We compared these methods of classification in terms of the accuracy percentage (= [number of samples misclassified/number of samples well-classified]*100).

Pathway analysis

To characterize the OPL subtypes, we performed an enrichment pathway analysis based on large-scale gene expression data in the four OPL datasets. Pathway-specific pathways were downloaded from the Molecular Signature Database (MSigDB database Molecular Signatures Database v5.2, 2016) and included a total of 1,329 canonical pathways. The ssGSEA was used to compute an enrichment score for each pathway in each sample[24,25] which was run from GenePattern.[26] Unlike GSEA which analyzes differential pathways between two phenotypical groups, the ssGSEA tool allows for computing an enrichment score (ES) of a given gene set in each sample. The gene expression values for a given sample are rank-normalized, and an ES is produced using the empirical cumulative distribution functions of the genes in the gene set and the remaining genes.[24,25] When UP and DN versions of a gene set are available, a combined score was computed. Default parameters were used (weighting component of 0.75 and minimal gene set size of 10 genes).

Deconvolution of the immune infiltrate in OPL

The Microenvironment Cell Populations-counter (MCP-counter) R package[27] was used for deconvolution of the immune infiltrate in each OPL sample. This method allows computing a score of the absolute abundance of eight immune cell types (T-cells; CD8 T-cells; Cytotoxic lymphocytes; NK cells; B lineage; Monocytic lineage; Myeloid dendritic cells; and Neutrophils) as well as of two stromal cell populations (endothelial cells and fibroblasts) in heterogeneous tissues from transcriptomic data. This score was defined as the log2 average expression of specific well-defined transcriptomic markers for each population.

Bioinformatics and statistics

Bioinformatics and statistics were performed using the Array Studio software (Omicsoft Corporation) as well as Bioconductor packages in the R language[56] and GraphPad Prism version 6.00 (San Diego, SA). Raw data from microarrays were processed using quantile normalization and the robust multi-array average (RMA) algorithm and were log2 transformed.[57] Custom CDF were used and downloaded from brainarray website.[58] miRNA raw data was log2 transformed for normalization. An integrative network including miRNA and differentially expressed genes between subtypes was built and visualized using miRNet[59] and the Gephi software,[60] respectively. The algorithm layout used in the Gephi software was ForceAtlas2. Hierarchical clustering analysis was done using the Pearson correlation and ward linkage method.

Unpaired Mann-Whitney or Wilcoxon tests were performed to compare the continuous values between the subtypes. Fisher’s exact test allowed us to analyze the contingency tables. In the 86 samples from the discovery set (GSE26549), oral cancer-free survival (OCFS) distributions were estimated using the Kaplan-Meier method and compared with the log- rank test. The OCFS was defined as the time from the first biopsy to oral cancer or to the date of last follow-up (for censored patients). The univariate and multivariate Cox proportional hazard model, including age, treatment arm (13cRA, RP with or without b-carotene BC) and histology (hyperplasia vs. dysplasia), were built to investigate the potential association between some biomarkers and OCFS in each molecular subtype. All statistical tests were two-sided, and P–values of 0.05 or less were considered to be statistically significant.