The FASTA file sequences for HPV16 and 18 E6/E7 proteins were obtained from UniProt (https://www.uniprot.org/). Online tools from Expasy (Swiss Institute of Bioinformatics), such as Protparam (https://web.expasy.org/cgi-bin/protparam/protparam), were used to analyze the molecular characteristics of these proteins, including molecular weight, theoretical isoelectric point, amino acid composition and average hydrophobicity. Protscale (https://web.expasy.org/protscale/) was used to predict the distribution of hydrophilic and hydrophobic regions in the proteins. Furthermore, the 3D structures of the HPV16 and 18 E6/E7 proteins were examined using the Phyre2 (http://www.sbg.bio.ic.ac.uk/~phyre2/) online software. These structures were modelled using Modeller (version 2.0; http://salilab.org/modeller/) and visualized using Rasmol (version 2.7.5.2; http://www.openrasmol.org/software/rasmol/) and Discovery Studio (2020; Dassault Systèmes SE; http://www.3ds.com/products/biovia/discovery-studio).

Raw gene expression data (datasets listed below) from HPV16 and HPV18-infected OSCC tissues were pre-processed to remove batch effects and normalize gene expression levels. Predefined gene sets from the GEO database and Molecular Signatures Database, were used to identify biological pathways and gene sets of interest. Genes were ranked based on their differential expression between HPV-positive and control (healthy) samples using metrics such as log-fold change or a moderated t-statistic. Gene set enrichment analysis (GSEA) was performed using software provided by the Broad Institute, Inc. (https://www.gsea-msigdb.org/gsea/index.jsp). Default settings were used unless otherwise specified. A focused search of the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) for oral cancer datasets using key words such as, ‘Oral cancer’, ‘HPV 16 and 18’ and ‘Homo sapiens’, and a filter for ‘Expression profile by microarray’ resulted in 11 potential datasets (21). In total, three datasets were selected after examining the requirement for complete clinical information, and excluding datasets containing data from blood samples, cancer cell lines or patients who had shown resistance to earlier therapies. Thus, three gene expression profiles [accession nos. GSE65858 (22), GSE42743 (23) and GSE6791 (24)] directly associated with OSCC and HPV were chosen from GEO datasets for further study.

Series matrix files for GSE65858, GSE42743 and GSE6791 were obtained for the analysis. Initially, the RNA-seq data from OSCC samples ensured that both HPV-positive and HPV-negative samples were obtained. The raw RNA-seq count data were normalized to obtain expression values using DESeq2 (Bioconductor) or edgeR (Bioconductor) (https://bioconductor.org/packages/release/bioc/html/edgeR.html). This file matched the order of samples in the gene expression matrix. Common nomenclature was used to convert the probe identifiers in each dataset into standardized gene symbols to ensure consistent gene identification (25). A robust multi-array average approach was used to normalize the datasets. Normalization was conducted utilizing R software (version 4.3.2; http://cran.r-project.org/bin/windows/base/), which standardized the gene expression data across all datasets by harmonizing them in scale and distribution (26). This process was crucial for mitigating systematic biases among samples or experimental conditions, facilitating accurate comparisons and analyses across datasets.

The GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) tool was used to identify DEGs within the OSCC and HPV datasets. The GEO2R tool produced a volcano plot that visually represents the variation in gene expression on the x-axis and statistical significance (P-value) on the y-axis, highlighting genes with significant changes in expression levels. DEGs were selected using the following criteria: P<0.01 and an absolute log-fold change >1. The FunRich (version 3.1.3; http://www.funrich.org/) tool was used to create Venn diagrams that showed the similarities and differences in DEGs across the three microarray datasets.

DEGs from OSCC samples were used to create a network that was explored using the STRING (https://string-db.org/) tool in the examination of PPIs within HPV-OSCC. Interactions with a confidence value of >0.4 were considered significant, which highlighted the reliability of these linkages (27). The network was visually depicted using Cytoscape (version 3.5.1; http://www.cytoscape.org), with the connections between proteins shown as lines of different thicknesses to illustrate interaction intensity. Hub genes were identified as proteins related to ≥10 others, which indicated their importance within the network. The present study utilized the MCODE (https://apps.cytoscape.org/apps/mcode) and cytohubba (https://apps.cytoscape.org/apps/cytohubba) plugins in Cytoscape to identify closely connected gene clusters. By setting specific parameters, such as a node score threshold of 0.2, a k-core value of 2 and a maximum depth of 100, the plugin could find and separate necessary modules or gene clusters that were closely linked.

Functional and pathway enrichment analysis is important for understanding the biological significance of identified DEGs and gene clusters (28). FOR GO analysis, the default enrichment statistic based on the Kolmogorov-Smirnov-like running sum statistic was employed. Gene set permutations were performed to assess statistical significance and calculate the false discovery rate (FDR). To ensure robustness, 1,000 permutations were used. Gene sets with FDR q-values <0.05 were considered significantly enriched. The Enrichr tool (https://maayanlab.cloud/Enrichr/) was employed to study biological functions and pathways associated with increased and decreased hub genes.

Additionally, a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed to identify pathways significantly enriched with the identified DEGs. To evaluate the significance of these pathways, criteria (normalization of DEGs) including a Benjamini-Hochberg adjusted P<0.05 were applied. A combined score based on the Jaccard coefficient (50%) and overlap coefficient (50%) >0.5 was considered statistically significant.

A total of 100 formalin-fixed paraffin-embedded (FFPE) tissue blocks were utilized, accessed from a tissue/sample databank held at the archives of Meenakshi Academy of Higher Education and Research University (Chennai, India). The tissue samples, which were collected over 5 years between May 2015 and June 2019, were obtained with ethics approval from the Institutional Review Board of Meenakshi Ammal Dental College and Hospital (approval no. MADC/IRB-XI/2017/235; February 7, 2023; Chennai, India). Patients had provided consent for their tissues to be used in future research. The sample set had 20 blocks from normal mucosa (NM) in healthy individuals (n=20), 40 blocks from patients with OL and 40 blocks from patients with OSCC. The age range of the healthy individuals was 32–50 years, with a mean age of 24 years and an age range of 18–29 years (male donors, n=10; female donors, n=10). Patients with OL had a mean age of 45 years, ranging from 32–50 years (male patients, n=28; female patients, n=12). Meanwhile, patients with OSCC were older, with a mean age of 50 years and an age range of 35–60 years (male patients, n=25; female patients, n=15). Gene expression levels of HPV16 and 18, along with their oncogenes E6 and E7, were detected using reverse transcription-quantitative PCR (RT-qPCR).

Sample preparation involved obtaining 10-µm slices from each FFPE tissue block. The sections were subsequently collected in 2 ml microcentrifuge tubes, ensuring an even distribution from the three groups. To remove the paraffin, 1 cc xylene preheated at 60°C for 10 min was added to each tube containing the tissue slices. The cells from the tissue pellet were lysed using a homogenizer, and the samples were centrifuged at a speed of 11,200 × g (deparaffinization at 60°C and inactivation at 95°C). During homogenization, Proteinase K (typically 20–50 µl, as per protocol) was added to digest the tissue. Once digestion was complete, samples were incubated at 95°C for 10 min to inactivate Proteinase K, ensuring no interference in subsequent analytical steps.

The DNA extraction process utilized a commercially available DNA isolation kit, QIAamp DNA FFPE Tissue Kit (Qiagen India Pvt. Ltd.). After the wax was removed, the tissue pellet was combined with 180 µl animal tissue lysis buffer and homogenized using a homogenizer. Subsequently, 20 µl proteinase K was added to the tubes containing the samples. The tubes were placed in an orbital shaking incubator at 56°C for 1–3 h until the tissue disintegrated. Next, 200 µl ethanol (96–100%) was added to the mixture. Following that, the mixture underwent 15 sec pulse-vortexing, then a brief centrifugation. The DNA extraction process utilized a commercially available DNA isolation kit, QIAamp DNA FFPE Tissue Kit (Qiagen India Pvt. Ltd.). After the paraffin wax was removed, the tissue pellet was combined with 180 µl animal tissue lysis buffer and homogenized using a homogenizer. Subsequently, 20 µl proteinase K was added to the tubes containing the samples. The tubes were placed in an orbital shaking incubator at 56°C for 1–3 h until the tissue fully disintegrated. Following this, 200 µl ethanol (96–100%) was added to the mixture, and the tubes were pulse-vortexed for 15 sec, followed by a brief centrifugation at 6,000 × g for 10 sec.

For RT-qPCR, cDNA was synthesized using the PrimeScript^™ RT Master Mix (Takara Bio Inc.) following the manufacturer's protocol, with incubation at 37°C for 15 min for cDNA synthesis and inactivation at 85°C for 5 sec. qPCR analyses were performed on the qPCR MX3000P system (Agilent Technologies Inc.) using the KAPA SYBR^® FAST qPCR Kit (KAPA Biosystems; Roche Life Science) with SYBR Green dye, which specifically binds to double-stranded DNA. Primers targeting HPV16, HPV18, 16-E6, 16-E7, 18-E6 and 18-E7 were used, with GAPDH as the reference gene, using forward (5′-TGCACCACCAACTGCTTAGC-3′) and reverse (5′-GGCATGGACTGTGGTCATGAG-3′) primers due to the stable expression of GAPDH across samples. Reactions were performed in triplicate with a no-template negative control. Thermocycling conditions included an initial denaturation at 95°C for 3 min, followed by 40 cycles of 95°C for 15 sec and 60°C for 30 sec. Melting curve analysis was conducted from 59 to 95°C to ensure product specificity. Quantification followed the 2^−ΔΔCq method (29), using Cq values for consistency. Samples were categorized as HPV-positive if Cq values for HPV16 or HPV18 primers fell below a set threshold, indicating detectable HPV DNA, and were confirmed by E6 and E7 region amplification, while samples with Cq values above the threshold were classified as HPV-negative.

Patients were categorized into HPV-positive and HPV-negative groups based on the presence or absence of HPV DNA, as determined by qPCR amplification of specific HPV genotypes (HPV16 and HPV18). If the Cq values for the HPV16 or HPV18 primers were below a predefined threshold, the sample was classified as HPV-positive. If no amplification was detected or if the Cq values were above the threshold, indicating no detectable HPV DNA, the sample was classified as HPV-negative. The HPV-positive cases were further confirmed by the amplification of the E6 and E7 regions of the HPV16 and HPV18 genomes, while the GAPDH gene served as a control for normalization.

Data concerning HPV16 and 18, and their E6 and E7 proteins, were consolidated into an Excel spreadsheet (version 16; Microsoft Corporation). Data were analyzed using SPSS (version 16; SPSS, Inc.). Fisher's exact tests were performed to ascertain statistical significance. P<0.05 was considered to indicate a statistically significant difference.

The physicochemical characteristics of HPV genotypes 16 and 18 E6/E7 oncoproteins were analyzed using the FASTA sequences sourced from UniProt (Table I). Hydrophilic and hydrophobic properties were assessed to understand their interactions in aqueous environments, which provided insights into their biological activities (Fig. 1A). Moreover, employing Modeller software, the study conducted a structural modeling of the E6 and E7 oncoproteins. The predicted 3D structures were achieved with a high degree of similarity, and these models were visualized using Rasmol, providing a clear depiction in Fig. 1B. This step offers a tangible representation of the proteins' potential structural conformations.

GSEA demonstrated upregulation of genes involved in critical pathways, such as the epithelial-mesenchymal transition (EMT), G₂/M checkpoint, inflammation regulation, apoptosis suppression, p53 and cancer-related pathways (Fig. 2A). Furthermore, Fig. 2B illustrates the binary relationship between HPV 16 and 18 E6/E7 oncoproteins, and their regulation of p53 and Rb proteins in human datasets. These findings strongly suggest that the primary focus of this study, which is the role of HPV 16 and 18 E6/E7 oncoproteins in the regulation of p53, contributes significantly to the process of oral carcinogenesis.

The DEGs between the OSCC and normal healthy tissues, as determined from the GEO datasets, are illustrated in Fig. 3A-C using volcano plots and mean difference plots, which indicate the differences in gene expression between the groups. The study suggests that HPV-associated OSCC may contribute to oral tumorigenesis through the regulation of specific pathways. Targeting these pathways could potentially yield improved therapeutic outcomes in the context of HPV-associated OSCC. This approach underscores the importance of understanding the molecular mechanisms underlying this type of oral cancer for the development of more effective treatment strategies.

DEG analysis across the datasets identified 1,145, 130 and 403 DEGs in the GSE6791, GSE42743 and GSE65858 datasets, respectively, all associated with HPV-related OSCC. This suggests potential as therapeutic targets or biomarkers, as modulating these genes may disrupt cancer-promoting pathways and enhance the efficacy of immunotherapies or other targeted treatments (Fig. 4A and B). The results of KEGG pathway analysis revealed significant enrichment in the ‘Interferon alpha/beta signaling’, ‘Interleukin-4 regulation of apoptosis’, ‘TSH regulation of gene expression’ and ‘EGFR1 pathway’ related to OSCC (Fig. 4C). Additionally, protein interaction predictions for HPV-associated OSCC were visualized using the STRING tool and a comprehensive network analysis was conducted with Cytoscape, utilizing plugins such as MCODE and cytoHubba (Fig. 4D). CytoHubba was used to identify hub genes from the PPI network, whereas MCODE was used to identify the genetic sequence for the specific proteins. This analysis unveiled the proteins associated with the p53 pathway from the OSCC datasets, highlighting their high nodal strength, closeness centrality, betweenness centrality and radiality, as detailed in Table II. These findings provide a comprehensive view of the molecular interactions and key players (EP300, HSP90AA1, TP53, CREBBP, NR3C1, SIRT1, MDM2, DNAJB1, H4C6 and STUB1) in the context of HPV-associated OSCC, which could be pivotal for further research and therapeutic development.

Heat map cluster analysis was conducted on DEGs from the GEO database and demonstrated that genes associated with HPV-related oral cancer were implicated in telomere regulation (Fig. 5A). GSEA plots (Fig. 5B) showed that HPV E6 and E7 genotypes may significantly impact the regulation, elongation and synthesis of telomeres, particularly on the lagging strand. Notably, these genes are associated with E6/E7 and play a role in influencing oral tumorigenesis, oncogenic pathways and telomeric regulation in oral cancer. These results provide compelling evidence of the connection between HPV-associated oral cancer and the regulation of telomeres, shedding light on the mechanisms underlying oral tumorigenesis.

The prevalence of HPV16 and 18, and their E6 and E7 oncoproteins, were compared across OSCC, OL and NM using Fisher's exact test (Table III). HPV16 E6 and E7 oncoproteins were significantly associated with OSCC, which indicated their higher prevalence in cancerous tissues. The HPV18 E7 oncoprotein was also significantly associated with OSCC. By contrast, no significant differences were found for the overall prevalence of HPV16 and HPV18. Additionally, the prevalence of HPV16 and 18, and their E6 and E7 oncoproteins, was compared among patients according to tumor location (Table IV), patient habits (Table V) and sex distribution (Table VI) among the OSCC, OL and NM groups. The analyses highlighted distinct patterns in the prevalence of HPV16 and 18, specifically in their oncoproteins E6 and E7, across OSCC, OL and NM groups. HPV16 oncoproteins E6 and E7 were significantly more prevalent in OSCC cases, with the E7 oncoprotein showing the highest association with OSCC, suggesting a strong link between HPV16 E7 and malignancy. Similarly, HPV18 E7 was significantly associated with OSCC, while no significant difference was found in the general prevalence of HPV16 or HPV18 across the three groups. In examining tumor location, a higher number of HPV-positive cases with E6 and E7 expression were found in OSCC at non-tongue sites compared with OL and NM. Regarding patient habits, HPV positivity (including E6 and E7 expression) was more prevalent among OSCC patients with tobacco use, notably in HPV16 E7-positive cases, further supporting an association between these oncoproteins and OSCC in patients with habits. Sex distribution data showed that males had a higher prevalence of HPV-positive cases across OSCC and OL groups, with no HPV positivity in NM, suggesting a sex-linked pattern in HPV infection prevalence within OSCC and OL.