The gene expression microarray datasets GSE3292, GSE6791 and GSE55542 were retrieved from the Gene Expression Omnibus (GEO) database. Only gene expression datasets that had been performed in HNC patients’ tumor samples and involved the HPV status were selected. In total 114 samples were analyzed, including 44 HPV⁺ and 70 HPV^- HNC patients ( Supplementary Table 1 ).

Quality control analysis, technical heterogeneity and batch effect exclusion, data background correction and data normalization are performed for all HNC samples. To detect the differentially expressed genes (DEGs) between HPV⁺ and HPV^- HNC patients, the LIMMA package of R was applied (27). Genes that met the cutoff criteria with fold change (FC) FC ≥ 2 and p-value (p) < 0.05 were considered as significant DEGs between HPV⁺ and HPV^- HNC patients. Overall patient’s transcriptome data, for a given gene, and the maximum expression value of the multiple probes were considered. The processed expression datasets were merged into a common expression matrix based on the gene names.

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using hypergeometric tests to uncover the putative biological roles of the DEGs between HPV⁺ and HPV^- HNC patients. Genes were annotated based on their associated Biological Processes (BP), Molecular Function (MF), and Cellular Component (CC) gene ontologies.

The Gene Set Enrichment Analysis (GSEA) was performed to identify gene sets that were differentially regulated between HPV⁺ and HPV^- HNC patients. Normal and absolute gene set enrichment analyses were conducted on the DEGs as stated before (28). The analyzed gene sets included the collections C2 (curated), C4 (computational), C5 (gene ontology), C6 (oncogenic signatures), C7 (immunologic signatures) and C8 (cell type signature) which are obtained on the latest updated version v2022 from the Human Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org). In addition to the subset of p53 target genes which was collected from online pathway databases, publications and through bioinformatics search ( Supplementary Table 2 ). The significant GSEA results were selected according to the nominal p < 0.05. Additionally, the leading-edge analysis was carried out to detect the biologically potential gene subset among the gene sets that are differentially regulated between HPV⁺ and HPV^- HNC patients.

To assess the potential Protein-Protein Interaction (PPI) network of the DEGs between HPV⁺ and HPV^- HNC patients, the STRING database (v11.5) was used (https://string-db.org/). The DEGs were mapped, and the corresponding PPI pairs were identified with the highest confidence interaction score of 0.9. Next, to identify the target genes, the PPI network was visualized using Cytoscape software (https://cytoscape.org/). In the network, nodes showing a higher number of links are apt to preserve the stability of the whole network based on the edge-connectivity concept. The cytoHubba application (v0.1) was used to find the top hub nodes ordered by degree of interactions using the Degree method.

To validate the protein expression of the selected potential genes from the in-silico transcriptomic analysis, a total of 65 patients with oral squamous cell carcinomas (OSCC) were collected, including 14 HPV⁺ and 51 HPV^- OSCC samples. The case classification of HPV status was based on the immunostaining of the p16 protein which is considered as a surrogate marker to diagnose high-risk HPV infection in OSCC according to the College of American Pathologists (29). All patients were selected and retrieved after ethical review and approval from the Tawam Hospital (Al-Ain, UAE) ethical committee (REC: AA/AJ/556). Detailed clinicopathological information was assembled from the patient’s medical records at Tawam Hospital. Two expert histopathologists (NA-R and AQ) initially examined the H&E-stained tissue sections to confirm the histological diagnosis. The tumors are originated from the tongue, floor of the mouth, cheek, gingiva, palate, or retromolar regions. Since the vermilion borders of the lip and the pharyngeal complex are not considered parts of the oral cavity, they were excluded.

The clinical and histopathological evaluation of the OSCC samples was performed based on the American Joint of cancer Committee on Cancer Staging (AJCC) classification (Sixth edition). All clinical information regarding the pathological tumor node metastasis (pTNM) staging of the original tumor, the level of cancer invasion and metastasis, the resection margin of cancer, the perineural invasion (PNI) and the lymphovascular invasion (LVI), were obtained from the patient’s histopathology reports.

Each tumor’s invasive front was identified, and a donor paraffin block’s core (0.5 cm) was considered for the experiment. From the 65 OSCC samples, 16-17 cores per block of tissue microarrays (TMA) were placed on four paraffin blocks. Initial H&E slides that subsequently revealed a limited tumor area were deleted from the following procedure analysis.

Paraffin sections were cut at 4 um thickness, placed on positively charged slides, and dried in an oven for 30 minutes at 70°C. Deparaffinization, rehydration, and target retrieval were performed in the PT Link (Dako) using a 3-in-1 procedure. Details of the immunohistochemistry (IHC) staining protocol for the used antibodies for p16, MDM2 and p53 were described in Supplementary Table 3 . IHC was performed on a Benchmark-ULTRA fully automated staining instrument (Roche Diagnostics, USA) using an Ultra View Universal DAB Detection kit from (Ventana). The antibodies were revealed by DAB and then counterstained with Meyer’s Hematoxylin and bluing reagent. Image acquisition and analysis were performed as previously described (30). Validation of Image J analysis for the three above-mentioned markers was performed by two observers who are experienced pathologists (NA-R and AQ). Controlling of threshold level was scaled and maintained without much adjustment (31). HPV positivity was considered when >40% of malignant epithelial cells were positive for P16 (32). The necrotic area of the specimen was not taken into account as previously described (31, 33).

The potential role of MDM2 in HPV-associated HNSC was further evaluated based on the Cancer Genome Atlas (TCGA) and the Genotype-Tissue Expression (GTEx) data resources using the Gene Expression Profiling and Interactive Analysis database (GEPIA) (http://gepia.cancer-pku.cn/), which represents the gene expression alterations in large cohort cancer samples, and the Tumor IMmune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/), which presents a comprehensive analysis of tumor-infiltrating immune cells and UALCAN resource (http://ualcan.path.uab.edu/index.html), which gives access to different publicly available cancer OMICS data. The survival analysis of the MDM2 gene was evaluated in HNSCC patients with high mutation burden, using the Kaplan-Meier plots.

The statistical analysis was performed using GraphPad prism software v5.0.1 (GraphPad Software, La Jolla California USA). Statistical significance between the protein expression and clinical parameters was analyzed using the t-test. Pearson’s correlation was applied to evaluate the possible correlation between protein expression and different clinical parameters across HPV⁺ and HPV^- OSCC patients. A two-way p < 0.05 was considered statistically significant.

After processing the whole transcriptomic data from the three preselected datasets GSE3291, GSE6791 and GSE55542, the merged data showed a total of 498 common differentially expressed genes (DEGs) (FC ≥ 2, p < 0.05) between the HPV⁺ compared to HPV^- HNSCC patients across the different datasets ( Figure 1A ). The DEGs included 207 down-regulated genes ( Figure 1B ) and 291 up-regulated genes ( Figure 1C ) in HPV⁺ compared to HPV^- HNSCC patients.

To assess the biological functions related to the identified DEGs, functional annotation and pathway analysis have been performed for GO and KEGG analyses. The up and down-regulated genes were separately annotated. The top 20 enriched terms including the related up/down-regulated genes in HPV⁺ compared to HPV^- HNSCC are presented in Figure 2 . The top 20 enriched terms and the related genes is provided in Supplementary Table 4 . The most up-regulated genes in HPV⁺ compared to HPV^- HNSCC were significantly enriched in the cell cycle regulation process including the PID E2F pathway (p= 1.89E-13), Transcriptional Regulation by E2F6 (p= 1.09E-09) and Regulation of cell cycle process (p= 9.99E-09), DNA replication and repair including DNA repair pathways (p= 2.06E-08), Chromatin organization (p= 3.6E-11) and Histone modification (p= 8.36E-08) and in the p53 tumor suppressor protein regulation (p= 5.02 E-05) including the Transcriptional Regulation by TP53 (p= 2.15E-09), PID P53 Downstream pathway and Regulation of cellular response to stress (p= 2.15E-09). While the down-regulated genes were mainly annotated in the VEGFA-VEGFR2 signaling pathway (p= 9.18E-09), Hemostasis (p= 6.71E-08) and Focal adhesion (p= 1.38E-07).

To identify the protein connections amongst the differentially expressed genes between HPV⁺ vs HPV^- HNSCC patients, the PPI-network interactions were constructed at the highest confidence of 0.9, defined as the minimum required interaction score ( Figure 3A ). Overall, the protein-protein interaction network has a p-value less than 1.0E-16, a total of 476 nodes with an average node degree of 0.929. The interaction enrichment showed a total of 221 edges while only 106 were expected, which indicates that the evaluated proteins are biologically interconnected as a group. For details of PPI networks and to identify the top protein complexes, we assumed that proteins under the same cluster probably share similar biological functions. Thus, the total proteins were clustered into the three highest-scored groups based on the clustering coefficient, using k-means clustering. Figure 3A showed that; Cluster 1 (colored in red), with PPI enrichment p-value = 1.13E-07, consisting of 27 nodes was highly enriched in focal adhesion and signaling by receptor tyrosine kinases. Cluster 2 (colored in green), with PPI enrichment p-value < 1.0E-16, consisting of 58 nodes was mainly associated with the DNA metabolic process and E2F regulation pathway. Cluster 3 (colored in blue), with PPI enrichment p-value = 2.97E-12, consisting of 38 nodes was related to cellular senescence and transcription regulation by TP53. The proteins and their node degree related to each cluster are listed in Supplementary Table 5 . Moreover, the PPI network was further explored to identify the hub genes, showing the highest number of connections in the network, using the Degree method under the CytoHubba plugin. The top 20 hub genes with high scoring were identified based on the degree of each protein node, including PCNA, MCM5, MCM3, MCM2, POLA1, RPA2, MCM6, PRIM1, RFC5, RBBP4, DUT, TYMS, RRM1, EZH2, LIG1, CHAF1A, RBBP7, MDM2, SUZ12, and SMC2 ( Supplementary Table 6 ). The sub-network interaction between the top 20 hub genes is presented in Figure 3B .

As anticipated, the GSEA showed that a subset of the p53-target genes was highly significantly enriched (ES = 0.558, p = 0.013 and FDR = 0.007) in the patients with HPV-related HNSCC ( Figure 4 ). The leading-edge analysis revealed that 18 core genes were significantly enriched in HPV-related HNSCC patients, including CDKN2A, DDB2, DUT, ERAP1, GLS2, MDM2, MLH1, NOTCH1, PCNA, POLH, PPM1D, PTTG1, SYK that were up-regulated in HPV⁺ HNSCC patients, however, CAV1, CD44, FLT1, ME1, and RTN4 were up-regulated in HPV^- HNSCC patients ( Figure 4 ).

To gain a better understanding of the potential distinctness in molecular signatures and biological mechanisms between HPV⁺ and HPV^- HNSCC. The total DEGs between HPV⁺ and HPV^- HNSCC patients were subjected to GSEA. First, the absolute GSEA was performed on 29,122 gene sets encompassing complex molecular mechanisms, different cellular pathways, biological processes and molecular functions. After excluding those with general terms or those including less than twenty genes, a total of 40 potent significant gene sets were differentially enriched between HPV⁺ and HPV^- HNSCC patients (p < 0.05). The details of these significant gene sets are listed in Supplementary Table 7 . Overall, in the identified significant gene sets, a broad overlap was noticed among the different gene sets with HPV-associated HNSCC-related cellular pathways or molecular mechanisms. Therefore, we grouped them mainly based on their involvement in the p53 pathway, Epithelial-Mesenchymal Transition, proliferation, GPR signaling and MAPK pathway, immune response/inflammation and metabolic process. Additionally, leading-edge analysis was carried out to detect the core-subset genes that potentially underlie the significant enrichment of the corresponding gene set and play substantial biological roles. Next, the gene frequency analysis was computed based on the occurrence of a gene across all the enriched leading-edge core genes from the 40 significantly over-represented gene sets. Based on the 50-percentile cutoff of the frequency values, 51 top frequent genes were identified ( Table 1 ). Interestingly, the p53-target genes were frequently listed in the core subset genes across several activated pathways, and it is noteworthy that the top two frequent genes across all significant gene sets are NOTCH1 (24 occurrences) and CAV1 (23 occurrences), which reinforces the previous analysis showing the important subset of p53-target genes in HPV-associated HNSCC.

To identify key biomarker genes and pathways involved in the progression of HPV-associated HNSCC, integrative analysis on the transcriptomics analyses was performed with several bioinformatics analyses. First, the comprehensive and integrated analysis was performed on the top significant pathways which resulted from the functional annotation of the DEGs, clustering of the top protein complexes and GSEA to show that the most marked observation to emerge from the data comparison was p53-related pathways and regulation. The second integrated analysis was carried out on the significant resulted genes from the top 20 hub genes analysis and the top frequent genes in GSEA to show intriguing correlation with five biomarkers which are EZH2, MDM2, PCNA, STAT5A and TYMS ( Figure 5A ). To further assess the expression stability of these genes across the different datasets, knowing that their ethnicities are different American, European-American and African American, the gene expression fold changes between HPV⁺ and HPV^- HNSCC were evaluated. Importantly, the MDM2 gene showed the highest gene expression difference between HPV⁺ and HPV^- HNSCC (Average log2FC = 1.89, Average of absolute deviation = 0.49) ( Figure 5B ). These results have further strengthened our assumption that over-expression of MDM2, proto-oncogene, may affect the occurrence and proliferation of HPV-associated HNSCC by disturbing the p53 target genes and consequently the p53-related pathways.

In order to provide additional evidence for our hypothesis, exhaustive in-silico validation of the MDM2 involvement in HPV-associated HNSCC was performed. First, the mRNA level of the MDM2 gene was assessed in a large sample size of HNSCC tumor samples (519 samples) compared to normal tissues (44 samples) using data from TCGA and GTEx in GEPIA. Importantly, the mRNA levels of MDM2 expression were significantly higher in HNSCC tissues than in normal tissues (p < 0.05) ( Supplementary Figure 1A ). In addition, the MDM2 expression was evaluated in HNSCC based on the HPV status using data from TCGA, including 41 HPV⁺ HNSC, 80 HPV^- HNSCC patients and 44 normal cases. As expected, the MDM2 showed a very significant over-expression in HPV⁺ than in HPV^- HNSCC patients (p = 2.39E-05) and in normal samples (p = 2.81E-07), however, only a slight difference was observed between HPV^- HNSCC patients and normal samples (p = 0.03) ( Supplementary Figure 1B ). These results were consistent with our findings from GEO datasets and further make evidence that the MDM2 gene is a key biomarker to distinguish between HPV⁺ and HPV^- HNSCC patients. Moreover, to assess the prognostic potential of MDM2, the overall survival analysis was performed in a cohort of HNSCC patients and in a cohort of HNSCC based on the HPV status using the Kaplan-Meier plots. The results reveal that a high level of MDM2 (HR= 1.44, 95%CI (1-2.06), p = 0.046) is correlated with poor overall survival in HNSCC patients ( Supplementary Figure 2A ). However, the high level of MDM2 demonstrates significantly increased overall survival (HR= 0.461, 95%CI (0.291-0.728), p = 0.002) in HNSCC HPV⁺ patients, while no significant association (p = 0.263) was reported between the MDM2 expression and the overall survival of the HNSCC HPV^- patients ( Supplementary Figure 2B ).

The clinical and histopathological evaluation of the collected OSCC samples was performed based on the AJCC classification (Sixth edition). Table 2 summarizes the characteristics of the 65 OSCC patients who were recruited in this study. The cohort included 70.8% males and 29.2% females. The patients ranged in age from 18 to 78 years old, with a median age of 57 and showing that about 61.5% of the patients were under the age of 60. Cigarette smoking habits constitute 21.5% of the patients while smokeless tobacco habits concern only 6% of the patients. The most common site of OSCC was the tongue (63.1%), followed by the cheek (13.8%) and Jawbones (9.2%). A total of 13.8% of the cohort had cancer in more than one place while about 64.6% of the cases were at the late stage (TNM stages III and IV). Recurrent tumors were observed in 21.5% of cases only. Lymph node metastasis was found in 35.4% of patients with single or multiple ipsilateral spread. About 50% of cases were well differentiated SCC with PNI in 28.6% of cases.

A tissue microarray was used for assessing the protein levels of MDM2 in HPV⁺ and HPV^- OSCC patients ( Figure 6A ). The IHC scoring showed a low level of IHC staining in HPV^- OSCC patients with an average expression of 37.87%, while, the HPV⁺ OSCC disclosed a high level of MDM2 expression with much more IHC staining indicating an average expression of 55.04% ( Figure 6B ). Importantly, according to the tissue microarray analysis, the MDM2 protein expression was significantly higher in HPV⁺ compared to HPV^- HNSCC patients (p = 0.031), ( Figure 6D ). Additionally, considering the close relationship between MDM2 and p53, we also evaluated the protein expression of p53 in HPV⁺ and HPV^- OSCC patients ( Figure 6C ). Notably, no significant variation (p = 0.44) has resulted between HPV⁺ and HPV^- OSCC patients ( Figure 6E ).

Mutual correlations were identified between MDM2 protein expression and different demographic clinical characteristics of the HNSCC patients, using Pearson’s correlation (p ≤ 0.05). As shown in Figure 7 , positive correlations have resulted importantly between MDM2 and HPV status, cardiovascular diseases, tumor size and tumor stage. Also, a slight positive correlation was reported between p53 and MDM2 expressions. While negative correlations have resulted between MDM2 and diabetes, smoking, metastasis and survival.