Unstimulated saliva was collected from clinically confirmed OSCC patients (n = 9) and normal control (NC) (n = 5) with no clinical evidence of OSCC. Saliva collection and RNA isolation was conducted as described previously (Pandit et al., 2013). Briefly, subjects were requested to refrain from drinking or eating for at least an hour before the collection of saliva. Approximately, 5 mL of saliva was collected from subjects into pre-chilled RNase-free conical tubes, and 5 μL of 20 U/μL RNase inhibitor (Ambion) was added immediately to prevent RNA degradation. Saliva was centrifuged at 2,600 g for 15 min at 4°C to remove cells and other debris, and the collected cell free saliva (CFS) was used for total RNA isolation by the TRIzol method (Ambion). Before saliva collection, informed consent was obtained from subjects, and all related protocols were approved by the Institute review board, Rajiv Gandhi Cancer Institute and Research Centre, New Delhi. The sample collection complied the Declaration of Helsinki.

RNA quantity and quality were analyzed on Bioanalyzer (Agilent 2,100). Small RNA-seq libraries were prepared using the Ion total RNA-seq v2.0 kit (Thermo Fisher Scientific) as per manufacturer’s protocol. To size-select the desired RNA products, Magnetic Bead Cleanup Module (Thermo Fisher Scientific) was used twice with the same sample. During the first bead binding, magnetic beads captured larger RNA species such as mRNA and lncRNA. During the second binding and with increased ethanol concentration, desired small RNA products (miRNA and other small RNA) in the supernatant re-bound to the magnetic beads. After washing the beads, the desired RNA products were eluted with pre-heated (80°C) Nuclease-free Water. Purified and size-selected RNA was converted into cDNA, and the barcodes and adaptors were ligated to the libraries. Library size distribution was verified by Bioanalyzer using a high-sensitivity DNA kit (Agilent) (Supplementary Figure S1). The libraries were diluted to 100 picomolar final concentration, and an equimolar pool was taken for clonal amplification. Template preparation for libraries on ionospheres was achieved using OneTouch 2 protocols (Thermo Fisher Scientific). The Ion 540 kit OT2 (Thermo Fisher Scientific) system in single-end mode with 200 bp chemistry was used for small RNA library sequencing.

Cutadapt (v1.8.1) algorithm was used to trim the 3′adaptor of the sequencing reads and assess the read quality. Sequencing reads were aligned to reference human genome (hg19) using Bowtie2 genome aligner (v2.2.5). Subsequently, mapped reads were screened for annotations using miRBase (v22) for miRNAs, and Gencode (v19) for other RNAs. Mapped sequencing reads were quantified using the Featurecount module (v1.4.6) of the R-package. For taxonomic classification of microbial reads, One Codex (Minot et al., 2015) algorithm was used. Total microbial diversity was calculated using the Shannon diversity index (read count ≥5, p-value ≤.05 was considered significant). The mean proportion of the sequencing reads aligned to bacterial taxa was calculated using Statistical analysis of taxonomic and functional profiles (STAMP) (v2.1.3) (Parks et al., 2014), p-value ≤.05. The differentially abundant taxa (p ≤ .05) were analyzed for associated metabolic pathways and gene ontologies (GO) using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.kegg.jp). The GOs associated with the metabolic functions were analyzed and visualized using Cytoscape (v3.8.2).

We used expression profiles of CFS RNAs from 5 NC and 9 OSCC samples. Since the gene expression matrix of the bodily fluids is often noisy, we selected only those genes that possess >5 read counts in at least 2 samples. Differential gene expression analysis was done using the DESeq2 web application (https://yanli.shinyapps.io/DEApp) Differentially expressed (DE) genes were further filtered based on their log₂ fold change (FC) ≥1 or ≤ -1, p-value <.05 for the downstream analysis. The volcano plots were built on the DESeq output using the GraphPad Prism 8.2.1.

To predict the functions of DE mRNAs and miRNAs target genes, we performed pathway enrichment analysis (cancer only) using the PANTHER (http://www.pantherdb.org) with a false discovery rate (FDR) < .05. The predicted pathways and pathway-related genes were analyzed and visualized as a network using Cytoscape (Shannon et al., 2003) (v3.8.2). Target genes of the differentially abundant miRNAs were predicted using Target Scan (v7.2) with a cumulative weighted context ++ score cut-off of -.4. The gene set enrichment analysis (GSEA) was conducted using iDEP.951 web application (http://bioinformatics.sdstate.edu/idep/) using Molecular Signatures (MSig) and oncogenic signatures reference gene sets. The GSEA was visualized as a bubble plot using SRplots (www.bioinformatics.com.cn/plot_basic_gopathway_bubbleplot).

The expression profile of DE genes was analyzed using Head and Neck Squamous cell carcinoma (HNSCC) dataset available on The Cancer Genome Atlas (TCGA) portal. The protein expression data was retrieved from HNSCC dataset available on Proteomics Data Commons (PDC) portal. Survival analysis was performed using Kaplan-Meier method on Km plotter web tool (https://kmplot.com/analysis/). The correlation of drug response with respect to DE genes was determined by Receiver Operating Characterization (ROC) method (https://www.rocplot.org). The DE miRNAs-mediated drug resistance pathways were predicted and visualized using ROC plotter (https://www.rocplot.org/mirna/index).

The statistical significance of DE genes was calculated using nominal p-value in DESeq2 web application. The significance in the difference between Shannon diversities of NC and OSCC was calculated using Student’s t-test. The enrichment in microbial taxa and metabolic pathways was calculated using STAMP statistical tool.

To identify the salivary RNAs that segregate OSCC patients from normal, we performed small RNA-seq of the saliva sampled from OSCC patients and NC (the clinicopathological characteristics are described in Supplementary File S1). Approximately, 27 % of reads from the NC and 20 % from OSCC aligned to hg19 (Alignment statistics are detailed in Supplementary File S2). Among the mapped reads, RNA biotypes including mRNAs, miRNAs, long intergenic non-coding RNAs (lincRNAs), small nuclear RNAs (snRNAs), transfer RNAs (tRNAs), and pseudogenes were detected in both OSCC and NC (Figure 1A). DESeq analysis of the sequencing data yielded a total of 166 RNAs that were DE between OSCC and NC, of which 145 were upregulated and 21 were downregulated, with a log₂ FC ≥ 1 or ≤ 1, p-value <.05 (Figure 1B; Supplementary File S3). Individual DE analysis of the RNA biotypes revealed 113 mRNAs (96 upregulated, and 17 downregulated) between OSCC and NC, (Figure 1C, top 30 DE mRNAs are shown in box). In addition, 7 miRNAs were found to be upregulated in OSCC compared to NC (Figure 1D) and 8 novel lincRNAs (5 upregulated, and 3 downregulated) (Figure 1E) are DE in OSCC. Notably, 35 novel tRNAs were found to be upregulated in OSCC (Figure 1F, top 10 tRNAs are shown in box). Other RNA species did not show any significant DE (Data not shown).

Salivary transcriptomes are enriched with oral microbial RNAs (Ostheim et al., 2020). Therefore, to examine whether the sequencing reads unmapped to hg19 were of microbial origin, we aligned these sequencing reads to the microbial genomes on ONE CODEX microbial analysis platform. Approximately, 9.76% of reads from OSCC, and 3.29% of reads from NC were aligned to annotated microbial genomes. Metagenomic analysis revealed a diverse range of bacterial phyla (Figure 1G) and genera (Figure 1H) in OSCC and NC samples. A large proportion of reads remain unmapped, we speculate that these reads plausibly belong to the other uncharacterized metagenomes whose sequences are unavailable to date. Summarily, we report an array of salivary RNAs that significantly segregate OSCC from NC, and RNAs annotated to oral microbiome.

The functional associations of DE mRNAs were examined by GSEA using gene sets of canonical oncogenic pathways available on MSig database. DE mRNAs were significantly (p < .05, FDR ≤.01) associated with multiple oncogenic signature pathways including BMl1, AKT, MTOR, KRAS, AKT, IL2, PGF and VEGF and p53 (Figure 2A, Supplementary File S4). In addition, GO analysis showed a significant association (p < .05) of mRNAs such as PRKCA, NOTCH2, ZAK3, CDH22, ERBB4, CSNK1E, PCDH15, CDHR3, ROCK1, and ADCY2 with signaling pathways that are known to be dysregulated in OSCC (Figure 2B). Since PRKCA was found to be associated with multiple oncogenic pathways, we further analysed its diagnostic and prognostic relevance using OSCC TCGA data. Notably, expression analysis revealed a significant upregulation of PRKCA at mRNA (Figure 2C) and protein levels (Figure 2D) in OSCC tumors compared to adjacent normal control (NAT). Survival analysis indicated that OSCC patients with higher expression of PRKCA were relatively correlated with low survival (Figure 2E). PRKCA activates signaling pathways involved in chemotherapeutic drug resistance (He et al., 2022), this agrees with our GO analysis, therefore, we examined the plausible role of PRKCA in therapeutic resistance of OSCC. Importantly, patients with elevated expression of PRKCA showed resistance to immune checkpoint inhibitors including anti-PD1, anti-PDL1 and anti-CTLA-4 for OSCC therapy, Figure 2F shows a combined immune checkpoint therapy response with respect to PRKCA expression in OSCC patients. Taken together, salivary PRKCA expression could be a potential diagnostic, prognostic, and therapeutic marker of OSCC.

The putative targets of DE-miRNAs were predicted using the Targetscan (Supplementary File S5) and subjected to functional enrichment analysis. The putative targets of miR-6087, miR-449b-5p, miR-3656, miR-326, miR-146b-5p and miR-497-5p were found to be significantly (p < .05) enriched in hallmark pathways of cancer (Figure 3A). Since the targets of DE miRNAs were associated with therapeutic resistance including p38 MAPK, EMT, RAS-PI3K, we investigated their putative role in chemo-resistance. Notably, patients with higher expression of miR-6087, miR-449-5p, miR-326 and miR-3656 showed no-response to platinum and/or taxane based drugs, compared to patients with lower expression of these miRNAs (Figures 3B–E). Conversely, the expression of miR-146b-5p and miR-497-5p was relatively lower in non-responders to platinum and taxane compounds (Figures 3F,G). Collectively, miR-6087, miR-449b-5p, miR-3656, miR-326, miR-146b-5p, and miR-497-5p can serve as diagnostic and drug resistance markers for OSCC.

Recently, oral microbial dysbiosis is strongly associated with the incidence and progression of OSCC (Lee et al., 2017). Since a significant proportion of our sequencing reads mapped to microbial genomes, we measured the differences in microbial composition between OSCC and NC. Shannon diversity index showed a significantly higher microbial diversity in OSCC compared to NC (p < .05) (Figure 4A). Differential enrichment of microbial types by STAMP analysis revealed that phyla Firmicutes, and Actinobacteria were significantly enriched in NC (p < .05), and Bacteroidetes were significantly enriched in OSCC (p < .05) (Figure 4B). Among the 15 genera, Prevotella was significantly enriched in OSCC (p < .05), and Neisseria and Haemophilus were enriched in the NC group (Figure 4C). GO analysis of differentially abundant microbial taxa revealed associations with pathways including inflammation, ROS generation, cell proliferation, amino acid metabolism, nucleic acid metabolism, and carbohydrate metabolism suggesting a plausible role in carcinogenesis (Figure 4D). Summarily, our investigation not only identified prominent microbial dysbiosis between OSCC and NC, but also highlights the importance of analysing microbial sequences in salivary RNA-seq studies.