Forty-eight FFPE tissue samples from patients diagnosed with OSCC, including both tumor and adjacent normal tissues from the southwest coast of India, were collected from Yenepoya Dental College and Hospital, Deralakatte within a duration of 2022–2024. The study received approval from the institutional ethical review boards of Yenepoya medical College and hospital ethics committee (Protocol No: YEC-1/2023/131). After quality control measures 13 samples were discarded. The remaining 35 samples were employed for the downstream analysis.

Whole exome sequencing (WES) was used to identify the mutation. Tumor and adjacent normal tissues of the OSCC patients were collected from the Yenepoya Dental College and Hospital for the mutational profiling. Out of the 35 samples 14 pairs (N = 28) were associated with the same patient (14 normal adjacent tissues as control and 14 OSCC tumor tissues as cases) from 14 individuals were employed in this study. Remaining seven were single tissue samples (five only tumor tissues and two only normal tissues). Total number of tumor tissues were 19 and total number of normal tissues were 16. The clinical details of the recruited samples is summarized in Table 1.

Whole Exome Sequencing (WES) libraries were prepared using Twist Exome 2.0 kit. Genomic DNA was extracted from the tissue samples, and samples with a DNA concentration of ≥20 ng/μL were selected for WES. All samples were sequenced on Illumina NovaSeq6000 system. The raw paired-end sequencing reads (PE150) generated from FFPE samples were initially processed for adapter removal and quality filtering to obtain high-quality trimmed reads. These processed reads were then mapped to the Homo sapiens reference genome GRCh37 using the Burrows–Wheeler Aligner (BWA) with default parameters (Li and Durbin, 2009). Variant calling followed standard best-practices workflows, where aligned BAM files were refined and variants were identified. The resulting VCF files were subsequently annotated using Ensembl’s Variant Effect Predictor (VEP) (Hunt et al., 2022) to assign functional consequences to each variant. A matched tumor–normal comparison strategy was used to distinguish somatic mutations from germline variants, emphasizing those alterations potentially driving tumorigenesis. All processed datasets comprised FASTQ, BAM, VCF, and annotated output files, as delivered in the sequencing report from NCGM.

To annotate the detected genetic variants, multiple clinically relevant databases and prediction tools were employed to classify variant types and assess their pathogenic potential. Tumor tissue was compared with adjacent normal tissue to distinguish somatic alterations from germline variants, particularly those implicated in tumor development. The dbSNP database was used as a comprehensive resource for cataloging genetic variations across different species and providing detailed information on each variant. ClinVar, a specialized repository for clinically important genetic alterations (Landrum et al., 2018), was utilized to identify variants classified as pathogenic or likely pathogenic. Additionally, integrated computational tools, PolyPhen-2 (Adzhubei et al., 2010) were applied to predict the functional consequences of the mutations. The OMIM® database was further consulted to determine any known links between the identified variants and genetic disorders.

To estimate PRS, we merged exome-wide data from 21 OSCC patients with genome-wide data of 204 South Asian individuals from our lab database using PLINK v1.9. The merged. bed file assessed 5,341 SNPs. PRSice2 (Choi et al., 2019) was employed to estimate PRS of oral cancer for the 21 OSCC patients and other individuals present in the merged dataset. Summary statistics from MVP_R4.1000G_AGR.Phe_145_2. META.GIA.dbGaP.txt was downloaded through GWAS catalogue and was employed as the reference (--base) data during PRS estimation. Notably, the reference file compared 2,203 oral cancer patients against 570,152 controls, assessing 17, 701, 051 SNPs. The test samples were divided into PRS quartiles; the lowest quartile (bottom 25%) was considered as reference. All subsequent PRS groups were considered to have an increased risk of developing oral cancer compared to the reference quartile.

This retrospective study analysed 48 patients with histopathologically confirmed Oral squamous cell carcinoma (OSCC) who underwent surgical resection at Yenepoya Medical College Hospital, Mangalore. Within a duration of (2022–2025) Formalin-fixed paraffin-embedded (FFPE) tissue samples were retrieved from the archives of Yenepoya Dental College Hospital after OSCC diagnosis confirmation. Collected clinicopathological data included (Table 2) age, gender, surgery date, primary tumor site, AJCC 8th edition TNM stage, tumor size, and tobacco use history.

The cohort comprised 18 (37.5%) females and 30 (62.5%) males. The mean age was 54.80 years (SD = 10.1), with a median of 55. The most common primary tumor site was the buccal mucosa (n = [10], (41.6%), followed by the [lateral border of tongue] (n = [6] % Choi et al., 2019) and the [lower alveolus and lower lip] (n = [5,3], [20.8,12.5] %)., TNM stage distribution, mean tumor size. The most common site was the buccal mucosa presenting as moderately differentiated OSCC (Table 3).

The table illustrates how tumor stage (pT1-pT4a) varies across different subsites in the oral cavity, highlighting the distribution of tumor advancement at diagnosis for each location. pT (Pathological Tumor): This refers to the size and extent of the primary tumor, determined after surgical removal and pathological examination. pT staging describes the size and local extent of the primary tumor after pathological examination. pT1-pT4a indicate increasing tumor size and/or invasion into nearby tissues, with specific measurements and involvement of adjacent structures defining each stage.

ClinVar classification system was used to determine clinically relevant genetic variants (SNVs and DELINs). Variants designated with classifiers such as pathogenic, likely pathogenic, risk allele, drug response, association, likely association, and protective were used for downstream analysis. We found discernible variation among the 19 recruited patients with tumor tissue samples in terms of the total number of clinically relevant mutations identified. While patient no.1 had the lowest number of (N = 132) clinically relevant variants, patient no.4 had the highest (N = 197).

Among the various types of genetic variants identified, the intronic variants comprised of 56% ± 7% of all identified variants, followed by missense variants (12% ± 1%) and non-coding transcript variants (10% ± 2%). Stop-gained mutations were found to be present in very low frequencies (<1%) (Figure 1).

The origin of the mutation was determined by comparing the mutations identified in the tumor tissue and those identified in the nearby normal tissue. Based on the origin of the mutations, out of an average of 171 clinically relevant mutations (range: 150–197), most were germline (93.3% ± 3.9%), and 6.7% ± 3.9% were somatic.

We identified 56 genes with at least one clinically relevant mutation in at least 80% of the recruited patients in the tumor tissue (Figure 2). Interestingly, mutations in 21 of these genes (ABCB1, CAPN10, CD44, CDC27, CHI3L1, CNN2, COMT, CTBP2, IGSF3, IL6, IL6-AS1, MAPKAPK3, MGST3, NAT2, NOS3, OPRM1, PADI2, SLC2A9, TBXT, TLR1 and VKORC1) were identified in all recruited patients. We identified 69 clinically relevant mutations that are present in at least 80% of the recruited patients in the tumor tissue (Figure 3). Among them, 19 mutations were identified among all recruited patients. The clinically relevant mutations that were identified in the tumor of at least 50% recruited patients are summarized in Supplementary Table S1.

We note here that, in The Cancer Genome Atlas (TCGA) - Head and Neck Squamous Cell Carcinoma (TCGA-HNSCC) dataset, recurrent somatic alterations are largely confined to well-established oncogenic drivers, including TP53, PIK3CA, and FAT1. In contrast, most genes found to be mutated in ≥50% of the recruited patients in our cohort predominantly encompassing drug-metabolizing enzymes, receptor genes, and immune-related regulators.

We identified exclusively tumor specific mutations by comparing the mutations identified in the tumor tissue with those identified in the nearby non-tumor normal tissue (N = 14 pairs). We identified ten mutations that were exclusively tumor specific, and present in at least 20% of the recruited patients (Table 4). Among them a single nucleotide variant rs5743618 of TLR1was found to be present in the tumor tissue of all recruited patients, followed by rs1805010 (intergenic) and rs546905091 (in PRKDC) that were present among 46% recruited patients. Notably, while the minor allele frequency (MAF) of rs5743618 is low (0.25) among European descents, it is very common among Asians (0.89 among South Asians and 0.99 among East Asians), underscoring the importance of multiple reference genomes to understand the global mutational landscapes (Karczewski et al., 2020). On the contrary, the pathogenic insertion mutation rs546905091 is extremely rare among South Asians (MAF = 0.0076), and this highlights its potential role in cancer (Auton et al., 2015). Interestingly, rs1805010, while designated as a pathogenic SNV in ClinVar, is also quite common among South Asians (MAF = 0.44) and also globally, indicating the importance of understanding the role of common genetic variants in the development of multifactorial disorders such as cancers.

Further, we systematically queried TCGA and focused on tumor-specific somatic variants, which revealed either very low or an absence of such mutations for these genes (Table 4). This disparity highlights a distinct mutational landscape in our cohort relative to TCGA, potentially reflecting cohort-specific biological features or differences in variant origin and detection.

Subsequently, a patient-wise comparison was performed to evaluate the overlap between mutations identified in this cohort and the most frequent canonical somatic drivers described in the TCGA-HNSCC dataset (Table 5). Of the established TCGA-HNSCC driver genes assessed, TP53 was the most frequently detected, present in five of 21 patients, followed by CDKN2A, identified in three patients. In addition, EGFR, a gene predominantly altered through copy-number amplification in TCGA-HNSCC, was observed in a single patient. Notably, several other genes that are commonly reported as recurrently mutated in TCGA-HNSCC: including FAT1, NOTCH1, KMT2D, NSD1, CASP8, and HRAS were not detected in this cohort. Collectively, these findings indicate a restricted overlap with the canonical TCGA-HNSCC mutational spectrum, suggesting that the genomic architecture of this cohort is shaped by distinct, potentially population- or exposure-specific factors rather than reflecting the dominant driver landscape reported in large, predominantly Western TCGA datasets.

Stratified analyses were performed to investigate whether the mutation landscape varied across key clinical and demographic dimensions, including age, sex, primary tumor site and pathological TNM stage (Supplementary Table S2). Comparison of mutation burden across disease stages using the Kruskal–Wallis test did not reveal statistically significant differences (p = 0.952), and post-hoc Dunn’s multiple comparison tests further confirmed that mutation burden did not vary significantly between any individual stage pairs (all adjusted p = 1.0). Likewise, comparison between node-negative and node-positive tumors (N0 vs. N+), used as a proxy for metastatic potential, demonstrated no significant difference in mutation burden (Wilcoxon rank-sum test p = 0.460). Subgroup analyses based on age category and sex also showed no discernible variation in mutational load across these strata. Fisher’s exact tests for each gene similarly did not identify any individual variant or gene disproportionately associated with advanced vs. early disease after false-discovery rate correction. Collectively, these results suggest that in this OSCC cohort, comprising patients with predominantly uniform gene panels, the global mutation burden and the presence of the pharmacogenomically associated variants assessed here do not appear to be strongly driven by TNM stage progression or major demographic parameters. Further, our results indicate that additional molecular mechanisms, such as copy-number changes, epigenetic dysregulation, immune microenvironmental factors, or environmental exposures, may have a more prominent role in OSCC clinical heterogeneity.

Out of the 40 most frequent genes identified in the tumor tissue that have at least one mutation in at least 95% of the recruited patients, 38 were found to be present in dbGENVOC, the comprehensive database with clinically pertinent somatic and germline variation data originating from >100 Indian oral cancer patients (Pradhan et al., 2021). This indicates that clinically relevant mutations identified in TBXT and PERM1, present in 100% and 95% recruited patients are novel mutations associated with OSCC patients from southwest India. Novel mutations in XRCC1 and ITPKB, previously identified by our team (Paul et al., 2025) was found among ∼80% and 95% of the recruited patients respectively, further highlighting the importance of understanding local mutational landscape. Notably, the missense variant rs25487 in XRCC1 has been associated with chemotherapy response. Among the 21 genes with exclusively tumor specific mutations in at least two recruited patients, 18 were found in dbGENVOC and two (CHROMR and DNAAF4) were novel. CHROMR (Cholesterol Induced Regulator of Metabolism RNA) (Barriocanal et al., 2022) is a long noncoding RNA (lncRNA) gene that is a key arbitrator of antiviral innate immune signaling in human while DNAAF4 (Dynein axonemal assembly factor 4) is involved in the preassembly of the multi-subunit dynein protein, which is essential for the proper functioning of cilia and flagella.

Gene ontology analysis was performed on the 21genes with exclusively tumor specific mutations in at least two recruited patients revealed. The biological processes (BP) of these genes are related to B cell lineage commitment, Replicative and Cellular senescence, Hematopoietic stem cell differentiation, Protein destabilization, rRNA processing, Cell cycle G1/S phase transition, and Regulation of G1/S transition of mitotic cell cycle. The main cellular components (CC) of these genes involve Striated muscle dense body, Cajal body, Nucleoid, Z-disc, I-band, Sarcomere, and Spliceosomal complex. The major molecular functions (MF) associated with these genes are Mitochondrial promoter sequence-specific DNA binding, Biotin transmembrane transporter activity, MDM2/MDM4 family protein binding, Cyclin-dependent protein serine/threonine kinase inhibitor activity, Disordered domain specific binding, Transcription coactivator binding, Protein kinase inhibitor activity and P53 binding.

The 21genes with exclusively tumor specific mutations in at least two recruited patients were analyzed for KEGG enrichment. The results revealed significant enrichment of these genes in various cancers and cancer-related pathways such as Pathways in cancer (hsa05200), Bladder cancer (hsa05219), Melanoma (hsa05218), Non-small cell lung cancer (hsa05223), p53 signaling pathway (hsa04115), Pancreatic cancer (hsa05212), Chronic myeloid leukemia (hsa05220) and Gastric cancer (hsa05226), as well as Platinum drug resistance (hsa01524) (Kanehisa et al., 2021; Gene Ontology Consortium, 2021) (Figure 4). In congruence with KEGG, Reactome pathway analysis revealed that these 21 genes are significantly associated with programmed cell death and senescence. This analysis also identified enrichment of these genes in immune system related pathways such as Interleukin signalling and disorders associated with immune system (Figure 5).

A protein–protein interaction (PPI) network of 21 genes with exclusively tumor specific mutations in at least two recruited patients was constructed to elucidate their interactions using STRING v12 (Figure 6). The network comprised of 19 nodes and 17 edges. The expected number of edges was 9. This indicates that our network has significantly more interactions than expected (p-value <0.00847), suggesting close interactions among these genes in various cancer related and TP53 associated biological pathways. Notably, AK2, NANOS1, DNAAF4, FRG1, IL4R, and MYO6 were not part of the PPI network.

The study participants, comprising OSCC patients and South Asian individuals from our lab database, were categorized into polygenic risk score (PRS) quartiles (Choi et al., 2019) using the lowest quartile (bottom 25%) as the reference group. All higher PRS quartiles were associated with an elevated risk of developing oral cancer relative to this reference (Q1 vs. Q2: Chi Square = 2.04, p-value = 0.15; Q1 vs. Q3: Chi Square = 4.15, p-value = 0.042; Q1 vs. Q4: Chi Square = 17.6, p-value = 2.72 x 10⁻⁵). The median PRS was −0.0055362, with an interquartile range (IQR) from −0.00595 to −0.0050734. Strikingly, 15 of the 21 OSCC patients were classified in the highest PRS quartile (Q4), suggesting a strong genetic predisposition to oral cancer (Figure 7). Of the remaining patients, four were in quartile 3 (Q3), and two were in quartile 2 (Q2). These findings underscore the potential of PRS as a valuable tool for cancer risk prediction.

A considerable subset of recurrently altered genes in OSCC are influenced by nuclear hormone receptor signaling, particularly pathways governed by the Vitamin D receptor (VDR), estrogen receptor (ER), and androgen receptor (AR). These receptors are critical regulators of epithelial cell differentiation, proliferation, and immune modulation, and their dysregulation has been linked to malignant transformation and tumor progression. To investigate potential mechanistic and therapeutic relevance, we evaluated whether the commonly mutated genes in our cohort are connected to these hormonal regulatory networks. The genes exhibiting consistent variation across all OSCC samples are summarized in Table 6, while the involvement of these genes in hormone receptor–mediated pathways is presented in Table 7.

We identified multiples genes that exhibited variations in all OSCC samples and have been documented to have associations with nuclear hormone receptor signaling (Table 8). Notably, several of these genes including ABCB1, CD44, CHI3L1, COMT, NAT2, NOS3, PADI2, TLR1, VKORC1, CYP19A1, CYP2D6, FKBP5, HMGCR, BGLAP, ESR1, CYP2B6, CYP3A5 and XBP1 — are linked to both Vitamin D receptor (VDR) and estrogen receptor (ER) pathways, highlighting the relevance of these hormonal axes in OSCC biology. Interestingly, all genes associated with VDR signaling also showed connections to ER regulation, and a subset further exhibited interactions with androgen receptor (AR) pathways. Given the established role of ER activation in promoting tumor progression and the antagonistic regulatory influence of VDR on ER signaling, these observations reinforce the importance of hormonal crosstalk in shaping oral carcinogenesis (Supplementary Table S3).

These findings further emphasize the need to investigate the functional implications of these variants in the context of the diverse environmental exposures and genetic backgrounds observed in the Indian population. Comprehensive characterization of how such hormone-associated molecular alterations influence tumor behavior may facilitate the discovery of actionable biomarkers and therapeutic targets. Our study supports the rationale for evaluating personalized treatment approaches that integrate vitamin D supplementation or hormone pathway modulators alongside conventional therapies to improve OSCC clinical outcomes.