This article is mentioned in:

Introduction

In recent years, there has been a significant increase in the incidence of oral malignancies associated with the human papillomavirus (HPV). Globally, HPV has been implicated in >38,000 cases of head and neck cancer, with ~30% classified as oropharyngeal malignancies (1). HPV types 16 and 18 are recognized as high-risk strains associated with the development of oral squamous cell carcinoma (OSCC). The infection status of HPV16 and HPV18 in OSCC shows that HPV16 is present in 20–35% of cases, while HPV18 is present in 5–10% of cases. Notably, India bears the highest burden of oral cancer in Asia, constituting ~20% of all reported cases in Asia, primarily due to delays in diagnosis that contributes to poor prognosis (2,3). HPV infections are estimated to contribute to ~20% of oral cancer cases and 60–80% of oropharyngeal malignancies, which often exhibit distinct clinical responses and survival rates. As such, HPV-positive tumors, especially oropharyngeal cancers, respond better to radiation and chemotherapy, likely due to increased radiosensitivity and enhanced immune response. Patients with HPV-positive tumors typically have higher survival rates, lower recurrence and reduced risk of metastasis compared with HPV-negative patients. This favorable prognosis has spurred research into less intensive treatment options for these patients compared with that for HPV-negative tumors (4). Besides HPV infection, several key factors contribute to OSCC. These include tobacco use (both smoking and smokeless forms), heavy alcohol consumption, chewing betel quid, poor oral hygiene, dietary factors (such as low intake of fruits and vegetables), sun exposure (particularly to the lips), genetic predispositions and immunosuppression, all of which increase the risk of developing OSCC (5). In addition to HPV16 and HPV18, other high-risk HPV types that can cause OSCC include HPV31, HPV33, HPV35, HPV39, HPV45, HPV51, HPV52, HPV56, HPV58 and HPV59. Although these types are less common compared with HPV16 and HPV18, they still contribute to the development of OSCC (6). The presence of these high-risk HPV types underscores the importance of regular screening and preventive measures, such as HPV vaccination, to reduce the incidence of HPV-related cancer (7). Understanding the broader spectrum of HPV types involved in OSCC can aid in developing more comprehensive strategies for prevention, early detection and treatment.

HPV is a non-enveloped virus known for its affinity for infecting epithelial cells. Its genome consists of two helical DNA strands within a spheroid structure, with the HPV family comprising ~100 unique types (8–10). Among the high-risk variants, HPV16 and HPV18 are most frequently implicated, in 50 and 20% of cervical cancer cases, respectively (11). Notably, these strains are responsible for a notable proportion of oral cancer cases among Indian betel quid chewers (12). The viral genome encompasses three distinct segments: Early (E), late (L) and upstream regulatory regions, with the E segment containing six proteins: E1, E2, E4, E5, E6 and E7, which aid in viral DNA replication and synthesis of new virus particles within infected cells (13). Of particular significance are the E6 and E7 oncoproteins, which serve a pivotal role in HPV-associated carcinogenesis by promoting cell proliferation, immortality and malignant transformation through interactions with key cellular proteins such as p53 and phosphorylated retinoblastoma protein (Rb) (14–16). Protein-protein interaction (PPI) networks and analysis of gene expression alterations associated with HPV infection provide insights into the molecular mechanisms underlying HPV-mediated oncogenesis (17–19). Structural analysis techniques can further demonstrate the alterations in E6 and E7 proteins as they interact with host proteins (20). The integration of data from various omics sources and employing bioinformatics tools can enable comprehensive exploration of molecular alterations implicated in HPV-associated cancer.

The primary aim of the present study was to validate the molecular profiles of HPV16 and 18 E6/E7 oncoproteins, and to elucidate their association with OSCC and oral leukoplakia (OL). This involved comparative analysis of differentially expressed genes (DEGs), functional annotations, pathway interrogation, and scrutiny of clinical and pathological indicators to detect high-risk HPV types, mainly focusing on the E6 and E7 oncogenes.

Materials and methods

Structural validation

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).

Data extraction

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.

Microarray data pre-processing

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.

Identification of differentially expressed genes (DEGs) in HPV-OSCC 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.

PPI network and module analysis of DEGs in HPV-OSCC 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.

Gene Ontology and pathway enrichment analysis

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.

In vitro studies

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).

Tissue sectioning and deparaffinization

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.

Extraction of DNA

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.

RT-quantitative PCR (qPCR)

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.

Statistical methods

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.

Results

Prediction of HPV genotypes 16 and 18 E6/E7 oncoprotein structures

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.

Figure 1. Structural validation of HPV genotypes 16 and 18 E6/E7 oncoproteins. (A) Secondary structure of HPV16 and 18 E6/E7 oncoproteins was demonstrated using the Protscale tool, which confirmed their hydrophilic and hydrophobic features. (B) HPV16 and 18 E6/E7 oncoproteins were examined in three dimensions using Rasmol and Discovery Studio software. HPV, human papillomavirus.

Table I. Structure prediction for HPV16 and 18 E6/E7 oncoproteins.

Table I.

Structure prediction for HPV16 and 18 E6/E7 oncoproteins.

HPV type	NCBI accession no.	Sequence	Length, aa	Molecular weight, kDa	Isoelectric point, pH	Aliphatic index
HPV16	NP 041325.1 (E6)	MHQKRTAMFQDPQERPRKLPQLCTELQ	158	19.18	9.16	68.48
		TTIHDIILECVYCKQQLLRREVYDFAFRDL
		CIVYRDGNPYAVCDKCLKFYSKISEYRH
		YCYSLYGTTLEQQYNKPLCDLLIRCINCQ
		KPLCPEEKQRHLDKKQRFHNIRGRWTGR
		CMSCCRSSRTRRETQL
	NP 041326.1 (E7)	MHGDTPTLHEYMLDLQPETTDLYCYEQL	98	11.02	4.20	78.57
		NDSEEEDEIDGPAGQAEPDRAHYNIVTFC
		CKCDSTLRLCVQSTHVDIRTLEDLLMGTL
		GIVCPICSQKP
HPV18	NP 040310.1 (E6)	MARFEDPTRRPYKLPDLCTELNTSLQDIE	158	18.87	8.95	78.99
		ITCVYCKTVLELTEVFEFAFKDLFVVYRDS
		IPHAACHKCIDFYSRIRELRHYSDSVYGDT
		LEKLTNTGLYNLLIRCLRCQKPLNPAEKLR
		HLNEKRRFHNIAGHYRGQCHSCCNRARQ
		ERLRRRETQV
	NP 040311.1 (E7)	MHGPKATLQDIVLHLEPQNEIPVDLLCHEQ	105	11.99	4.70	71.59
		LSDSEEENDEIDGVNHQHLPARRAEPQRH
		TMLCMCCKCEARIELVVESSADDLRAFQQ
		LFLNTLSFVCPWCASQQ

[i] HPV, human papillomavirus.

GSEA of HPV-associated pathways

GSEA demonstrated upregulation of genes involved in critical pathways, such as the epithelial-mesenchymal transition (EMT), G2/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.

Figure 2. Gene Set Enrichment Analysis detection of HPV-related OSCC. (A) Gene expression in OSCC tissue samples compared with normal tissue across the datasets GSE65858, GSE42743 and GSE6791, depicting increased transcript levels in OSCC. The color gradient represents the transition from cancerous to healthy samples. (B) Interaction of HPV16 and 18 E6 and E7 oncoproteins and their roles in OSCC through analysis on the UniProt server, demonstrating their contribution to OSCC. Blue indicates the presence of p53 in HPV 16 and 18 (E6 and E7). HPV, human papillomavirus; OSCC, oral squamous cell carcinoma.

DEGs in OSCC

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.

Figure 3. Analysis of DEGs in HPV-associated OSCC. DEGs in the OSCC Gene Expression Omnibus datasets (A) GSE65858, (B) GSE42743 and (C) GSE6791, are shown using volcano plots (top) and the mean-difference approach (bottom). Upregulated genes, red; downregulated genes, blue; non-significant DEGs, black; HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; DEGs, differentially expressed genes.

Identification of key pathways in HPV-associated OSCC

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.

Figure 4. Identification of genes associated with HPV-induced OSCC. (A) FunRich software was used to create a Venn diagram, showing 12 common genes across OSCC samples from the GSE65858, GSE42743 and GSE6791 datasets. Enriched (B) biological processes and (C) KEGG pathways linked to HPV-OSCC. (D) Illustration of protein interactions in these GSE datasets using Cytoscape. (E) Network density using MCODE and (F) identification of key genes with strong interactions using Cytohubba in the OSCC protein network. HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table II. Targets in human papillomavirus-associated oral squamous cell carcinoma datasets.

Table II.

Targets in human papillomavirus-associated oral squamous cell carcinoma datasets.

Top DEGs	Betweenness	Closeness	Clustering co-efficient	Degree	Radiality	MCODE	MCC method
EP300	5.20	11.0	0.45	22	2.27	7.00	8,640
HSP90AA1	5.20	11.0	0.43	22	2.27	7.00	8,640
TP53	5.20	11.0	0.43	22	2.27	7.00	8,640
CREBBP	3.41	10.5	0.43	20	2.18	6.37	7,200
NR3C1	1.80	10.0	0.41	18	2.09	6.37	7,200
SIRT1	1.78	10.0	0.41	18	2.09	6.37	7,200
MDM2	1.23	10.0	0.40	18	2.09	6.37	6,480
DNAJB1	1.23	10.0	0.40	18	2.09	6.00	6,480
H4C6	0.33	09.0	0.37	14	1.90	5.78	1,440
STUB1	0.28	09.0	0.36	14	1.90	5.78	1,440

[i] Genes were a part of high-risk HPV sub type. DEG, differentially expressed gene; MCC, Maximal Clique Centrality; MCODE, Molecular Complex Detection.

E6/E7 expression profile associated with cellular signaling

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.

Figure 5. Key cancer-causing and telomeric signaling pathways in HPV and oral cancer datasets. (A) Heat map comparing OSCC and normal samples demonstrating the substantial rise in oncogenic and telomeric pathway regulator-interacting gene expression. Blue denotes the low intensity and red indicates the high intensity of tumor in the samples. (B) Gene Set Enrichment Analysis plots demonstrating the effects of HPV-associated OSCC on cancer-causing sites in telomere extension, polymer switching, telomere extension and telomeric C strain lagging strand synthesis. Enrichment was measured using the NES. HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; NES, net enrichment score.

In vitro analysis

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.

Table III. Prevalence of HPV16 and 18, and E6 and E7 oncoproteins, in OSCC (n=40), OL (n=40) and NM (n=20) groups.

Table III.

Prevalence of HPV16 and 18, and E6 and E7 oncoproteins, in OSCC (n=40), OL (n=40) and NM (n=20) groups.

HPV type and oncoprotein	OSCC, n (%)	OL, n (%)	NM, n (%)	P-valuea
HPV-16	3 (7.5)	5 (12.5)	0 (0)	>0.999
E6	7 (17.5)	2 (5.0)	0 (0)	0.0152
E7	8 (20)	0 (0)	0 (0)	0.0000229
HPV-18	6 (15)	6 (15)	0 (0)	0.1577
E6	0 (0)	0 (0)	0 (0)	-
E7	6 (15)	0 (0)	0 (0)	0.0070

a Fisher's exact test. HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; OL, oral leukoplakia; NM, normal mucosa.

Table IV. Positivity of HPV16 and 18, and E6 and E7 oncoproteins, according to location on the tongue and other sites among OSCC (n=40), OL (n=40) and NM (n=20) groups.

Table IV.

Positivity of HPV16 and 18, and E6 and E7 oncoproteins, according to location on the tongue and other sites among OSCC (n=40), OL (n=40) and NM (n=20) groups.

A, Tongue
		HPV16			HPV18
Group	No. of cases	HPV16 positive	E6	E7	HPV18 positive	E6	E7
OSCC	13	2	2	3	0	0	1
OL	8	1	0	0	3	0	0
NM	0	0	0	0	0	0	0
B, Other sites
		HPV16			HPV18
Group	No. of cases	HPV16 positive	E6	E7	HPV18 positive	E6	E7
OSCC	27	1	5	5	6	0	5
OL	32	4	2	0	3	0	0
NM	20	0	0	0	0	0	0

[i] HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; OL, oral leukoplakia; NM, normal mucosa.

Table V. Positivity for HPV16 and 18, and E6 and E7 oncoproteins, in patients with or without tobacco habits among OSCC (n=40), OL (n=40) and NM (n=20) groups.

Table V.

Positivity for HPV16 and 18, and E6 and E7 oncoproteins, in patients with or without tobacco habits among OSCC (n=40), OL (n=40) and NM (n=20) groups.

A, With habits
		HPV16			HPV18
Group	No. of cases	HPV16 positive	E6	E7	HPV18 positive	E6	E7
OSCC	25	1	3	7	2	0	3
OL	34	4	2	0	4	0	0
NM	0	0	0	0	0	0	0
B, Without habits
		HPV16			HPV18
Group	No. of cases	HPV16 positive	E6	E7	HPV18 positive	E6	E7
OSCC	15	2	4	1	4	0	3
OL	6	1	0	0	2	0	0
NM	20	0	0	0	0	0	0

[i] HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; OL, oral leukoplakia; NM, normal mucosa.

Table VI. Distribution of HPV-positive and negative cases in OSCC (n=40), OL (n=40) and NM (n=20) groups according to sex.

Table VI.

Distribution of HPV-positive and negative cases in OSCC (n=40), OL (n=40) and NM (n=20) groups according to sex.

A, HPV positive
Group	No. of cases	Male	Female
OSCC	17	12	5
OL	10	9	1
NM	0	0	0
B, HPV negative
Group	No. of cases	Male	Female
OSCC	23	16	7
OL	30	29	1
NM	20	10	10

[i] HPV, human papillomavirus; OSCC, oral squamous cell carcinoma; OL, oral leukoplakia; NM, normal mucosa.

Discussion

HPV, particularly HPV16, has emerged as a key etiological factor in a subset of OSCC cases. Unlike traditional risk factors, such as tobacco and alcohol use, HPV-associated OSCC predominantly affects younger individuals, often with no history of smoking or heavy drinking (30). Understanding the molecular mechanisms underlying HPV-driven OSCC is important for early detection, prognosis and targeted therapies. The primary aim of the present study was to investigate the structural features of HPV16 and 18 E6/E7 oncoproteins through a physicochemical investigation. Modeler software was used to successfully create 3D structural models demonstrating how the E6 and E7 oncoproteins were spatially arranged. The identification of DEGs in HPV-associated OSCC, including EP300, HSP90AA1, TP53, CREBBP, NR3C1, SIRT1, MDM2, DNAJB1, H4C6 and STUB1, demonstrated critical insights into the molecular mechanisms of the disease. These genes may serve pivotal roles across various pathways essential for cancer development, such as the EMT, G2/M cell cycle checkpoint, cancer-associated pathways, regulation of inflammatory responses, suppression of apoptosis and the p53 pathway. EMT is important in cancer metastasis, as it enables epithelial cells to gain migratory and invasive properties (31). Through epigenetic modifications, genes such as EP300 and SIRT1 regulate transcription factors pivotal for EMT (32), which illustrates the complexity of cancer metastasis at the epigenetic level. The G2/M checkpoint of the cell cycle ensures DNA integrity before mitosis, with TP53 serving a critical role in its activation in response to DNA damage (33). The balance between TP53 and its negative regulator, MDM2, is key for cell cycle control (34), pointing to potential areas of intervention in cancer treatment. HSP90AA1 helps in protein folding by synthesizing cellular protein folds, thereby helping the tumor cells evade cellular stress and increasing tumor survival (35). CREBBP and EP300 modify transcription factors needed for growth and proliferation (36). The improper translation of these two genes can cause cells to grow and survive without control, which are fundamental characteristics of cancer cells. NR3C1 is involved in anti-inflammatory reactions (37), whereby controlling inflammation could be a therapeutic target. Apoptosis suppression, facilitated by the interplay between TP53 and MDM2, and the role of SIRT1, underscores the evasion of programmed cell death as a cancer hallmark. The p53 pathway, central to DNA damage response and cellular stress, is disrupted in HPV-OSCC, particularly by the HPV E6 protein promoting p53 degradation (38). This disruption highlights the significance of the p53 pathway in maintaining cellular homeostasis and preventing cancer development. This suggests that the degradation of p53 by HPV E6 in HPV-OSCC disrupts the cellular stress response and impairs the ability of cells to repair DNA or undergo apoptosis, promoting oncogenesis and therapeutic resistance.

The results from the KEGG pathway analysis indicated significant enrichment in the ‘Interferon alpha/beta signaling’, ‘Interleukin-4 regulation of apoptosis’, ‘TSH regulation of gene expression’ and ‘EGFR1 pathway’ in the context of HPV-associated OSCC, thus highlighting the intricate network of cellular signaling cascades implicated in the pathogenesis of this disease. The IFN-α/β signaling pathway has an important role in the innate immune response to viral infections, including HPV (39). IFNs are cytokines that induce antiviral states in cells, upregulating genes that inhibit viral replication and spread. TSH-regulating gene expression serves an important role in metabolic regulation and cellular proliferation (40), which indicates that its involvement in OSCC reflects the complex hormonal and metabolic reprogramming that occurs during cancer progression. The EGFR1 signaling pathway is instrumental in cell growth, differentiation and survival (41). The significant enrichment of the EGFR1 signaling pathway in HPV-OSCC highlights the importance of growth factor signaling in the development and progression of this type of cancer. HPV oncoproteins are known to interact with and dysregulate growth factor pathways to promote cellular proliferation and evade apoptosis (42). PPI network analysis further elucidated the extensive interactions of proteins within HPV-linked OSCC datasets. This comprehensive examination of the molecular interaction's sheds light on the key nodes and factors influencing the development of oral cancer. Telomere regulation analysis (through in vitro pathway study with molecular databases analysis) demonstrated the potential role of E6 and E7 oncoproteins in stabilizing telomeres, which suggests an association with unregulated cell division in oral cancer. The results of in vitro assays that assessed HPV presence across diverse tissue types reinforced the epidemiological association between HPV infection and the onset of oral malignancies. This examination corroborated the molecular findings from other study analyses, providing a more holistic depiction of HPV involvement in OSCC (41,43). Fundamentally, the present findings substantiated the association between HPV infection and OSCC, emphasizing the significance of recognizing HPV as a pivotal determinant in the prevention, diagnosis and management of oral malignancies. Recognition of the inherent constraints, such as potential biases in datasets, underscores the imperative for ongoing validation and expansion of these findings.

The present study identified unique gene expression profiles and molecular pathways in HPV16 and HPV18-infected OSCC tissues, which differed significantly from non-HPV-related OSCC. HPV-positive samples exhibited upregulation of E6 and E7 oncogenes, which may disrupt the expression of tumor suppressor genes, such as p53 and Rb, due to viral protein interactions. By contrast, in non-HPV-related OSCC cases, these genes are often mutated or altered through other carcinogenic pathways (43). The molecular pathways activated in HPV-positive OSCC are distinct, with the PI3K/AKT pathway and cell cycle regulation prominently affected by viral oncoproteins (44–47). By contrast, non-HPV-related OSCC shows alterations in pathways associated with tobacco and alcohol exposure, such as the EGFR pathway and oxidative stress responses.

Clinically, patients with HPV-positive OSCC generally have an improved prognosis and respond differently to treatment modalities, such as radiation and chemotherapy, compared with those with non-HPV-related OSCC (48,49). These differences underscore the importance of distinguishing between HPV-related and non-HPV-related OSCC in clinical management and treatment planning. Furthermore, the present research identified specific biomarkers and potential therapeutic targets unique to HPV16 and HPV18 infections. Therapeutic strategies targeting the HPV oncoproteins E6 and E7 could provide more effective treatments for HPV-positive OSCC, while different approaches might be necessary for non-HPV-related cases.

Through structural validation, gene expression analysis and in vitro experiments, the present study demonstrated how HPV, particularly types 16 and 18 and their E6 and E7 oncoproteins, contribute to oral cancer development. These findings highlight the dysregulation of key pathways, such as EMT, cell-cycle control and inflammation, shedding light on the underlying mechanisms of HPV-associated OSCC. These results provide important implications for targeted therapies and diagnostic strategies in HPV-related oral cancer, underscoring the need for further research to advance the current understanding and clinical management of these diseases. Future studies should incorporate protein interaction network analyses utilizing platforms such as STRING and Cytoscape to elucidate the complex molecular interactions and pathways implicated in these processes. Such analyses will provide both visual and quantitative insights into the underlying molecular dynamics. Moreover, it is crucial to investigate the specific roles and functions of these key proteins within the identified networks. Conducting thorough experimental validation will strengthen the current findings and may also reveal new therapeutic targets. Such future studies are essential for translating this research into clinical applications for oral cancer, aiming to enhance treatment strategies, improve patient outcomes, and achieve more effective and lasting tumor suppression and tissue regeneration. In conclusion, the increasing incidence of mouth cancer linked to HPV, especially in the oropharyngeal area, is a substantial global health concern. The main aims of the present study were to investigate the structural properties of HPV16 and 18 E6/E7 oncoproteins, and to identify associations between HPV genotypes 16 and 18 and both OSCC and OL. The analysis verified the presence of high-risk HPV types 16 and 18, together with their corresponding E6/E7 oncogenes, in both OSCC and OL through RT-qPCR analysis on DNA isolated from FFPE tissue samples. The findings demonstrated increased occurrence of HPV18 in both the OSCC and OL groups in comparison to HPV16. Moreover, further investigation is required to identify small inhibitors that may be potential targets for treating oral cancer linked with HPV.

Acknowledgements

Not applicable.

Funding

The authors would like to thank the Deanship of Scientific Research at King Khalid University for supporting this work through Large Group Project (no. RGP-2/504/45).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

SFJH, SSA, KB, LT, IM, IHM, NSA and MIK contributed to the conception and design of the present study. Data collection was performed by SFJH, KB, LT and IM. Data analysis was performed by SSA, LT. and IHM. The first draft of the manuscript was written by SFJH, KB, LT, IM and MIK. Reviewing and editing of the manuscript was performed by SSA, IHM and MIK. SFJH and IHM confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

The present study was conducted in accordance with The Declaration of Helsinki and was approved by the Institutional Review Board of Meenakshi Ammal Dental College and Hospital (Chennai, India; approval no. MADC/IRB-XI/2017/235; February 7, 2023).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References