Identification of keratinocyte-associated genes for immune characterization and drug response prediction in oral squamous cell carcinoma

Data sources

The Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) dataset was accessed to obtain the scRNA-seq data, clinical information, and survival information of 12 OSCC samples (GSE215403). In addition, the dataset of OSCC microarray data (GSE30784) was also downloaded and screened to obtain bulk RNA-seq data of 45 normal oral tissue samples and 167 OSCC tissue samples.

Processing and analysis of the scRNA-seq data of OSCC

Quality control and cell filtering of the scRNA-seq data of OSCC were conducted using the “Seurat” package (Butler et al., 2018). A total of 31,776 cells were retained after eliminating those with 200 to 6,000 expressed genes and more than 10% mitochondrial content. After normalizing and scaling the raw data using the SCTransform function, principal component analysis (PCA) was performed (Bahudian & Valdovinos, 2025). The “harmony” package (Korsunsky et al., 2019) was employed to remove the batch effect amongst various samples, followed by conducting uniform manifold approximation and projection (UMAP) for dimensionality reduction using the RunUMAP function. Next, we built a K nearest neighbors (KNN) network based on the Euclidean distance with the top 50 principal components employing the FindNeighbors function. The cells were then clustered into subpopulations by the FindCluster function at the resolution of 0.05. The intercellular communication of keratinocyte subpopulations mediated by ligand–receptor pairs was analyzed using the R package “CellChat”(Efremova et al., 2020), with “Cell-Cell Contact” as the interaction type. Finally, the results were visualized into bubble plots.

Analysis of specific high-expressed genes and functional annotation between cellular subpopulations

The FindAllMarkers tool (parameters: “logfc.threshold = 0.30, min.pct = 0.30, only.pos = T”) was applied to screen specifically high-expressed genes (significance threshold of p < 0.05) from various subpopulations as cellular markers. Functional annotation of differentially expressed genes (DEGs) was performed using the package “CusterProfiler” (Yu et al., 2012).

High-dimensional weighted gene co-expression network analysis screening of keratinocyte-related genes

According to the characteristics in scRNA-seq data, high-dimensional weighted gene co-expression network analysis (hdWGCNA) can be employed to explore particular gene expression patterns involved in various biological processes (Langfelder & Horvath, 2008) and build co-expression networks across several cellular and spatial hierarchical scales (Morabito et al., 2023). Here, we used the hdWGCNA to read single-cell transcriptome rds data to develop the co-expression network of keratinocytes under the optimal soft threshold of eight. The modules related to keratinocytes were sectioned by calculating the connectivity of the modules, and the key genes within each module were selected based on their connectivity to the modules for further study.

Identification of DEGs

DEGs between the normal and OSCC tissue samples in the GSE30784 dataset were filtered using the “limma” package (Ritchie et al., 2015) under the criteria of |log₂FC| ≥ 1 and padj < 0.01. Then, genes present in intersection between the selected modular genes and the DEGs were considered as differentially expressed keratinocyte-related hub genes.

Construction of a gene diagnostic model

We first tested the diagnostic potential of the hub genes using the R package “multipleROC” (Cook, 2008). Subsequently, the R package “rms” was used to construct a diagnostic model for predicting the risk of OSCC, and a nomogram was drawn (Iasonos et al., 2008). The discriminative performance of the model was compared based on the receiver operating characteristic (ROC) and area under the ROC curve (AUC). The accuracy and efficiency of the model were validated by plotting calibration curves with the “caret” package, and decision curve analysis (DCA) (Gerds, Andersen & Kattan, 2014) was conducted utilizing the “rmda” R package.

Immune infiltration analysis

The official CIBERSORT website (https://cibersortx.stanford.edu/) was accessed to source the expression data of 22 types of common immune infiltrating cells (LM22), and the level of immune cell infiltration between OSCC samples and control samples in the GSE30784 dataset was calculated using the “CIBERSORT” R package (Newman et al., 2015). Rank-sum test was applied to determine significant immune infiltration differences (p < 0.05) between the two types of samples, while Spearman’s correlation coefficient (p < 0.05) was used to calculate the relationship between immune cell infiltration and crucial genes. Using the R language’s “estimate” package (Yoshihara et al., 2013), immune infiltration in the GSE30784 dataset was examined and the results were displayed as corresponding scores. The correlation between the pivotal genes and ESTIMATEScore, StromalScore, and ImmuneScore was analyzed according to the Spearman’s correlation coefficient, and the rank-sum test was employed to determine whether the immune scores were significantly different between the two types of samples (p < 0.05).

Cell culture and plasmid transfection

Human normal oral keratinocytes (HOK; CP2610) were ordered from ScienCell Research Laboratories (Carlsbad, CA, USA). Human OSCC cells (HSC-3; JCRB0623) were purchased from the Japanese Collection of Research Bioresources (JCRB, Osaka, Japan). Human oral squamous carcinoma cells (WSU-HN30; STM-CL-5612) were purchased from Stemmer Biotechnology Co (Shanghai, China). DMEM medium (Gibco, Thermo Fisher Scientific Inc., Waltham, MA, USA) supplemented with 1% penicillin-streptomycin (15140122; Gibco) and 10% fetal bovine serum (FBS; 10099141C, Gibco) was utilized for cell culture in an incubator with 5% CO₂ at 37 °C. The cells were identified by Short Tandem Repeat (STR) analysis, and the results of mycoplasma detection for these cells were negative.

According to the instruction of Lipo3000 Liposome Transfection Reagent (L3000-001; Thermo Fisher Scientific Inc.), WSU-HN30 and HSC-3 cells (2 × 10⁴ cells/well) in the logarithmic growth phase were transfected with KRT16 knockdown plasmid (si-KRT16, sequence: 5’AACAGCGAACTGGTACAGAGC 3’) and control plasmid (si-NC), which were purchased from GenePharma (Shanghai, China).

RNA extraction and quantitative real-time polymerase chain reaction

Total RNA was extracted from HOK, WSU-HN30, and HSC-3 cells using the RNA Extraction Kit (TRIzol, 15596026; Invitrogen, Carlsbad, CA, USA), following the protocols. The purity and concentration of the total RNA were determined, and cDNA template was generated using the HiScript II First-strand cDNA Synthesis kit (R211-01; Vazyme, Nanjing, China) (Zhang et al., 2023). The quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using specific primers and the KAPA SYBR^® FAST kit (KK4600; Sigma Aldrich, Burlington, MA, USA). The data were calculated using the 2^−ΔΔCT method with GAPDH as the internal control. See Table 1 for the primer sequences of the specific genes.

Cell viability

According to the protocol, CCK-8 (CK04; Dojindo, Tokyo, Japan) was performed to assess the effect of KRT16 on the viability of WSU-HN30 and HSC-3 cells. Briefly, the cells (5 × 10³ cells/well) were grown in 96-well microtiter plates for 24, 48, and 72 hours (h). After washing the cells twice with phosphate-buffered saline (PBS), 100 µL of fresh medium and CCK-8 solution (10 µL) was added to each well for 3-h incubation at 37 °C with 5% CO₂. Absorbance (at 450 nm) was measured in a SPECTROstar^® Nano (BMG LABTECH GmbH, Ortenberg, Germany) (Ma et al., 2023).

Cell migration assay

For wound healing assay, the transfected cells were inoculated into 6-well plates (5 × 10⁵ cells/well) and 2 ml of cell suspension was added for incubation in an incubator at 37 °C with 5% CO₂. When the cells adhered to the wall, the monolayer was wounded with a 10 µL plastic pipette tip to form a uniform scratch. PBS was used for washing the monolayers, which were then incubated in medium without FBS. The wound edge distances between the two edges of the migrating cell sheet were photographed at 0 h and 48 h. All the experiments were conducted in triplicate.

Cell invasion assay

WSU-HN30 and HSC-3 cell suspensions (5 × 10⁵ cells/well) were prepared in serum-free medium. Then, 100 µL of cell suspension was supplemented to the upper Transwell chamber (Corning, Corning, NY, USA) pre-coated with Martigel (30 µg/well; BD, San Jose, CA, USA), while the lower chamber contained 600 µL of DMEM medium and 10% FBS. The migrated cells were fixed by 4% paraformaldehyde and colored by crystal violet solution. The cells in the lower chamber were counted under a microscope from six different fields of view.

Cell apoptosis assay

For the flow cytometry assay, the cells were rinsed in phosphate buffered saline and suspended using 1× binding buffer provided by the assay kit (C1062S, Beyotime, China) to adjust to the concentration of 1 × 10⁶ cells/mL. 100 µL cell suspension was then taken and added with the working solution of Annexin V-FITC and propidium iodide (5 µL) for a 10-minute mixture reaction in the dark at ambient temperature. The apoptosis was finally tested in a flow cytometer (BD, USA).

Statistical analysis

GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA) and R software version 3.6.0 (R Foundation, Vienna, Austria) were utilized in all statistical analyses. Data were presented as mean ± standard deviation and analyzed using the Student t-test or two-way ANOVA followed by Bonferroni correction for post hoc multiple comparisons. A p < 0.05 stood for statistically significant difference.

Identification of cell subtypes in OSCC

We obtained a total of 12 cell clusters after processing the scRNA-seq data of OSCC using the “Seurat” package (Fig. 1A), and eight major cell subpopulations were identified, namely, T cells, keratinocytes, macrophages, B cells, fibroblasts, plasma cells, endothelial cells, and mast cells (Fig. 1B). Subsequently, cell type annotation was performed using the CellMarker2.0 database to identify representative markers for T cells (CCL5, GNLY, and GZMA), keratinocytes (KRT17, KRT14, and S100A8), macrophages (LYZ, C1QB, and C1QA), B cells (CD37, CD79A, and MS4A1), fibroblasts (DCN and LUM), plasma cells (IGHG1 and IGHG4), endothelial cells (VWF and PECAM1), and mast cells (TPSB2 and CPA3). The normalized expression levels of these marker genes were displayed as bubble plots (Fig. 1C). Keratinization of the oral mucosal epithelium has been reported to affect the immune response, and immune checkpoint-associated factors are low-expressed in OSCC of keratinized epithelial origin (Kitsukawa et al., 2024). To further explore the potential impact of keratinocytes on OSCC, we performed enrichment analysis on specifically high-expressed genes of the keratinocyte subpopulation. It was found that these genes were significantly enriched in establishment of protein localization to endoplasmic reticulum entries, translational initiation, signal-recognition particle (SRP)-dependent cotranslational protein targeting to membrane, and endoplasmic reticulum-targeted protein (Fig. 1D, p < 0.05).

Cellular communication between keratinocytes and other cells in OSCC

The reciprocal receptor–ligand pair communication network between keratinocytes and other cell subtypes in OSCC was explored employing CellChat analysis, which showed that keratinocytes communicated with other cellular subtypes. The ligand–receptor information for each cluster was extracted, and we discovered that endothelial cells affected keratinocytes through SELE-CD44 and APP-CD74, fibroblasts affected keratinocytes through APP-CD74, and macrophages affected keratinocytes through HLA-DRB5-CD4 (Fig. 2). These results indicated that keratinocytes and cell subpopulations had a high level of reciprocal communication.

Keratinocyte-related modules screened by hdWGCNA

To develop a co-expression network, a soft threshold power of β = 8 was found to achieve an optimal network connectivity when the scale-free topology fit reached 0.90 (Fig. 3A). Additionally, nine modules (Keratinocyte-M1 to Keratinocyte-M9) were identified (Fig. 3B) by calculating the expression of module eigengenes and module connectivity. The top 10 core genes of each module were selected based on the connectivity of the characterized genes (Fig. 3C).

Subsequently, we calculated the expression of genes in each module in different cells. Notably, the genes in the M3 and M5 modules showed high levels of activation mainly in keratinocytes (Fig. 4A), therefore these two modules were regarded as the key modules. Then, correlation analysis that M6 was negatively correlated with the M3 and M5, while strong positive correlations were observed between the M3, M4, M7, and M9 as well as between the M1, M6, and M8 (Fig. 4B). Then, co-expression network for the two modules were mapped using their top 25 hub genes (Figs. 4C–4D), with the 10 most highly connected genes in the inner circle as key genes and 15 genes in the outer circle as secondary key genes.

DEG analysis and screening of signature genes between OSCC and control samples

The “limma” package screened 1573 DEGs (789 upregulated genes and 784 downregulated genes) between OSCC and control samples (Fig. 5A). Next, the intersection of the hub genes of the two key modules and the DEGs contained four genes (KRT6B, KRT16, CSTB, and CSTA), which were considered as differentially expressed keratinocyte-related genes (Fig. 5B). The expression heatmap of the four genes showed that KRT6B and KRT16 were high-expressed in OSCC samples, while CSTB and CSTA genes were low-expressed (Fig. 5C). A diagnostic model for OSCC was developed using the four genes, and the model reliability was assessed according to the ROC curves and AUC. It was found that the AUC of KRT6B, KRT16, CSTB, and CSTA reached an AUC of 0.791, 0.917, 0.95, and 0.902, respectively (Fig. 5D), with KRT16, CSTA, and CSTB showing an AUC greater than 0.9. This indicated a strong diagnostic value of the three genes in OSCC, and these genes were employed in subsequent analysis.

Diagnostic model construction

To further improve the risk prediction for OSCC patients, we build a nomogram combining the expression of the three genes (Fig. 6A). The model showed that CSTB had the greatest impact on OSCC prognosis, followed by KRT16, while CSTA had limited impact on the prediction of the OS. Next, the ROC curve showed that the AUC of the model was 0.996 (Fig. 6B), implying that the diagnostic model was highly reliable. The calibration curve demonstrated that the calibration curve was close to the standard one, which suggested that the nomogram had strong prediction performance (Fig. 6C). In addition, the reliability of the model was tested utilizing the DCA decision curve. It can be observed that the gain of the nomogram was notably higher than the gain of other individual genes, suggesting a strong predictive ability of our model (Fig. 6D).

Correlation of the three characterized genes with immune infiltration

Immune cell infiltration differences between OSCC samples and control samples in the GSE30784 dataset were compared. The results of CIBERSORT showed that the infiltration of most immune cell types was notably different between the two types of samples. In particular, cells such as macrophages M1 (M1), macrophages M0 (M0), and mast cells activated had higher levels of infiltration in OSCC patients than in control samples, while B cells memory, T cells CD8, immune cells such as macrophages M2 (M2), and mast cells resting had lower levels of infiltration (Fig. 7A). Subsequently, we calculated the correlation between the expression of the characterized genes and the immune cell infiltration scores. As shown in Fig. 7B, the expressions of the characterized genes were linked to the infiltration of a majority of immune cells. Specifically, the expression of CSTA and CSTB was positively linked to the infiltration scores of B cells memory, T cells CD8, M2, dendritic cells resting, and mast cells resting (p < 0.01) but negatively correlated (p < 0.001) with the infiltration of cells such as M0, M1, and mast cells activated. In contrast, the expression of KRT16 was positively associated with the infiltration scores of mast cells activated, M0, and M1 (p < 0.01) but negatively associated with the infiltration scores of B cells memory and T cells CD8 (p < 0.001). Immune microenvironment differences between OSCC samples and control samples were also compared. ESTIMATE analysis revealed that StromalScore, ImmuneScore, and ESTIMATEScore were all significantly higher in OSCC samples than in control samples (p < 0.001), indicating a higher level of immune cell infiltration in the OSCC samples (Fig. 7C). We also found that the expression of CSTA and CSTB was negatively correlated with StromalScore, ImmuneScore, and ESTIMATEScore in OSCC (p < 0.001), while KRT16 expression was positively linked to ESTIMATEScore and StromalScore (Fig. 7D, p < 0.05). These results indicated that the characterized genes may be closely associated with immune cell infiltration in OSCC.

Downregulation of KRT16 suppressed the migratory and invasive abilities of OSCC cells

The mRNA expressions of three key genes in HOK, WSU-HN30 and HSC-3 cells were detected by performing qRT-PCR. The results demonstrated that the expression of KRT16 was significantly upregulated in WSU-HN30 and HSC-3 cells than in HOK cells, while the expression of CSTA and CSTB was significantly lower in WSU-HN30 and HSC-3 cells than in HOK cells (Fig. 8A). As we have confirmed a high diagnostic value of KRT16 and the potential role of its overexpression in the malignant progression of OSCC, we prioritized KRT16 for functional validation. After the knockdown of KRT16 gene in OSCC cells (Fig. 8B), it was observed that KRT16 knockdown significantly inhibited the survival of WSU-HN30 and HSC-3 cells (Fig. 8C). Subsequently, wound healing and transwell assays showed that KRT16 knockdown markedly suppressed the migration and metastasis of WSU-HN30 and HSC-3 cells (Figs. 8D–8G) but notably increased apoptosis in WSU-HN30 and HSC-3 cells (Fig. 8H). These results indicated a potential role of KRT16 in the occurrence of OSCC and its development.