Datasets containing samples of both cSCC and normal skin tissue were downloaded from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). Four datasets based on the GPL571 platform, including GSE42677 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42677), GSE45164 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE45164), GSE66359 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66359), and GSE98767 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98767), were finally selected for this study. After removing samples, in which systemic treatments, such as immunosuppression or chemotherapy, were used, a total of 46 cSCC samples and 27 healthy skin tissue samples were included in the subsequent analysis (Table 1). DEGs between the cSCC and control groups were identified in the downloaded data matrices using the R package LIMMA (version 3.40.6; R Software for Statistical Computing, Vienna, Austria). The screening criteria were set as follows: fold change (FC) $>$2 and p-value 0.05. The results were visualized using volcano plots and heat maps. The Venn diagram was drawn using an online Draw Venn Diagram tool (https://bioinformatics.psb.ugent.be/webtools/Venn/).

Gene Ontology (GO) functional enrichment analysis, pathway enrichment analysis (including the Kyoto Encyclopedia of Genes and Genomes, Reactome Gene Sets, and WikiPathways), and disease enrichment analysis based on the DisGeNET database were used, and data were plotted as bar graphs using METASCAPE (https://metascape.org/gp/). A threshold of the enrichment number of genes $\ge$3 and a p-value 0.01 were selected.

Interactions between the DEG-encoded proteins were analytically mapped using STRING 12.0 (https://cn.string-db.org/) with a threshold of a combined score $>$0.4 (i.e., medium confidence level). The obtained protein-protein interaction network was analyzed using MCODE, a CYTOSCAPE plugin (https://apps.cytoscape.org/apps/mcode) for the acquisition of functionally relevant key modules with a degree cut-off $>$2, node score cut-off $>$0.2, k-score of 2, and maximum depth of 100. The resulting module genes were verified by another plug-in, CYTOHUBBA (https://apps.cytoscape.org/apps/cytohubba), utilizing the Matthews Correlation Coefficient (MCC) algorithm.

The SIGnaling Network Open Resource (SIGNOR3.0), a repository of manually annotated causal relationships between human proteins, biologically relevant chemicals, stimuli, and phenotypes, provides advanced graphical tools to explore human cellular networks. By searching the hub genes in the key modules using the SIGNOR3.0 database (https://signor.uniroma2.it/), we created networks of potential signaling pathways. The causal relationships and mechanisms of how these proteins in the pathway act on the development of cSCC were further determined by searching the literature and analysing the molecular structure of functional domains obtained from the cBio Cancer Genomics Portal (cBioPortal).

The expression levels of TPX2, KIF11, CDK1, and pituitary tumor transforming gene (PTTG1) in squamous cell carcinoma of skin and other organs were validated using Gene Expression Profile Interaction Analysis (GEPIA; http://gepia.cancer-pku.cn/) based on The Cancer Genome Atlas (TCGA) dataset. Kaplan–Meier plots were generated online (http://kmplot.com/analysis/index.php?p=service&cancer=pancancer_rnaseq) to plot overall survival curves based on the expression of four core genes in head and neck squamous cell carcinoma (HNSC); the best cutoff value was selected automatically. The Human Protein Atlas (HPA) (https://www.proteinatlas.org/) provided immunohistochemical expression maps of protein expression patterns in cSCC and normal skin tissues.

Genetic information, such as mutation rate, prognostic information, molecular structure and mutations, and post-translational modification sites for patients with cSCC were obtained from cBioPortal (https://www.cbioportal.org/). Samples were downloaded from two studies: Cutaneous Squamous Cell Carcinoma (MD Anderson, Clin Cancer Res 2014) and Cutaneous Squamous Cell Carcinoma (UCSF, NPJ Genome Med 2021).

Small molecule compound-gene interaction analyses were based on the Comparative Toxicogenomics Database (CTD) (downloaded in November 2016), and the network was visualized using the NetworkAnalyst website (https://www.networkanalyst.ca/).

The SDF format of NSC689534 was downloaded from the PubChem database, subsequently imported into ChemDraw 3D for energy minimization using the MM2 module to obtain the lowest energy advantage concept, and then saved as a mol2 file. Protein structures were downloaded from the UniProt database and visualized separately using the PyMOL Molecular Graphics System (Version 2.5, Schrödinger, LLC). Afterward, the ligands (NSC689534) and receptors (TPX2 and KIF11) were saved as pdbqt files using MGTools (Version 1.5.6, Center for Computational Structural Biology, Scripps Research, USA) after processing for dehydration, hydrogen addition, charge calculation, merging of non-polar hydrogens, and others. The ligands and receptors were docked using AutoDock Vina 1.1.2, [17, 18] with the docking parameters shown in Table 2. Then, concepts with higher scores were visualized using PyMOL and Discovery Studio (BIOVIA, Dassault Systèmes, FR). Pharmacokinetics, drug-likeness, and medicinal chemical analyses of NSC689534 were performed using the SwissADME algorithm. We analyzed drug-likeness on the basis of the Lipinski (Pfizer Inc., New York, NY, USA) rule of five [19], whereas the Abbot bioavailability score and gastrointestinal absorption were used to estimate the oral bioavailability of the drug [20].

The A431 (National Collection of Authenticated Cell Cultures, Beijing, China) and HaCaT (Newgainbio, Wuxi, China) cells were grown in Dulbecco’s modified Eagle’s medium (SH30022.01, HyClone, Logan, UT, USA) supplemented with 10% fetal bovine serum (SH30074.03, HyClone) in 5% CO${}_2$ at 37 °C. Cell lines were authenticated with Short Tandem Repeat profiling (MicroreaderTM21 ID System, Applied Biosystems®, Waltham, MA, USA). Mycoplasma test was performed in all cell lines every other week using the Mycoplasma Real-Time PCR Detection Kit (Cat. No.4460623, MycoSEQTM Applied Biosystems®, Warrington, Cheshire, UK).

Cells were washed twice in cold phosphate buffered saline and, then, lysed in radioimmunoprecipitation assay buffer for 30 min. After being boiled for 10 min, equal amounts of protein were separated via SDS-PAGE and transferred to polyvinylidene fluoride membranes. The membranes were blocked with 5% skim milk and 0.1% Tween 20 (P1379, Sigma-Aldrich, St Louis, MO, USA) in tris-buffered saline with Tween-20 (TBST) and were then incubated overnight at 4 °C with antibodies against KIF11 (23333-1-AP, Proteintech, Rosemond, IL, USA) and TPX2 (11741-1-AP, Proteintech) in 5% skim milk. After being washed, the membranes were incubated in secondary antibodies for 2 h. Protein bands were detected using an enhanced chemiluminescence system.

The anti-cancer ability of NSC689534 was assessed using A431 cell line. 15,000 cells/well were seeded in 96-well plates for 24 h, treated by relevant concentration (0, 0.25, 0.5, 0.75, 1, 2.5, 5, 7.5 µM) of NSC689534 (obtained from MedChemExpress, Monmouth Junction, NJ, USA) for 8 h or 24 h. The cells were then incubated with fresh complete media containing 10% cell counting kit-8 (CCK-8) reagent (C6005, NCM, Suzhou, CHN) for 2 h. The absorbance at 450 nm was measured with the microplate reader manufactured by Bio-Tek Instruments (Winooski, VT, USA). Every group was repeated 3 times.

RNAs were extracted from HaCaT and A431 cells with 1 mL TRIzol (Cat. No.10296010CN, TRIzol™ LS reagent, Carlsbad, CA, USA) per sample. cDNA synthesis and qPCR were performed using HiScript® Ⅲ RT SuperMix for qPCR (R323, Vazyme, Nanjing, Jiangsu, China) and AceQ Universal SYBR qPCR Master Mix (Q511, Vazyme). GAPDH was tested for sample loading control. Primer (No.2514129018, Sangon, Shanghai, China) sequences were as follows: TPX2: TPX2-F GTCACCAAATCAGTTGACTTCC, TPX2-R AAGGATGCTTTCGTAGTTCAGA; KIF11: KIF11-F CATACTCTAGTCGTTCCCACTC, KIF11-R CAACCAAGTTCAACTTTCCGAT; GADPH: GAPDH-F CAGGAGGCATTGCTGATGAT, GAPDH-R GAAGGCTGGGGCTCATTT.

Analyses and graphs were conducted with GraphPad Prism 9.4.1 (GraphPad Software, Inc., San Diego, CA, USA). One-way ANOVA was used for data analysis in western blotting while T-test in RT-qPCR. p 0.05 was statistically significant.

We downloaded conforming microarray datasets from the NCBI GEO database (Table 1) and examined shared genes that were differentially expressed between cSCC and normal skin tissue in all datasets. All DEGs identified are shown as volcano plots (Fig. 1a–d). Based on the matrices of the GSE42677, GSE45164, GSE66359, and GSE98767 datasets, 3082, 1000, 1332, and 12327 DEGs, respectively, meeting the threshold of $\mid$FC$\mid$ $>$2 and p 0.05 were screened. Of these, 51 DEGs were shared among the four datasets (Fig. 1e).

Fig. 1.

Identification of 51 differentially expressed genes (DEGs) from Four Datasets. (a–d) A total of 51 DEGs identified from four datasets. Differentially expressed genes between cutaneous squamous cell carcinoma (cSCC) and healthy tissue samples in the datasets GSE42677 (a), GSE45164 (b), GSE66359 (c), and GSE98767 (d). DEGs were screened based on FC $>$2 and a corrected p-value 0.05. Red data points in the volcano plot represent upregulated genes whereas blue points represent downregulated genes. (e) Venn diagram of common DEGs among the four datasets. cSCC, cutaneous squamous cell carcinoma; DEG, differentially expressed gene; FC, fold change.

The 51 DEGs were significantly enriched in squamous cell carcinoma of the skin, actinic keratosis, and other skin-related disorders (Fig. 2a). GO and Reactome Gene Set analysis showed enrichment in cornified envelope/formation of cornified envelope (Fig. 2b,c), indicating that the screened DEGs were closely related to skin diseases. Pathway enrichment analysis also demonstrated that these 51 genes were associated with a variety of squamous cell carcinoma-related signaling pathways, including human papillomavirus infection, transforming growth factor beta signaling pathway, tyrosine kinase receptor signaling pathway, and regulation of the TP53-related pathway (Fig. 2c).

Fig. 2.

Enrichment analysis of the 51 DEGs. (a) Bar chart of significant GO terms with DEGs. (b) Bar chart of significant enriched pathways. (c) Bar chart of significant enriched diseases. The bar graphs were plotted using METASCAPE with a threshold of the enrichment number of genes $\ge$3 and a p-value 0.01. DEG, differentially expressed gene; GO, gene ontology.

The PPI network comprizing these 51 genes contained 28 nodes and 53 interactions (Fig. 3a), and their GO functional enrichment analysis was mainly related to peptide cross-linking, keratinization, epidermal cell differentiation, and skin development. Using the MCODE plug-in of CYTOSCAPE, two functionally relevant key modules were extracted from the PPI network (Fig. 3c,d). Module A contained eight genes and 23 edges, whereas module B contained five genes and 10 edges. Nearly all 13 hub genes were highly expressed in the four GEO datasets (Fig. 3b). Genes in Module A were enriched in keratinization, peptide cross-linking, and keratin-forming cell differentiation, while those in Module B were mostly associated with spindle organization and polarization, mitotic spindles, and the cell cycle. In addition, we found that the scores of the genes in both modules were ranked in the top 15 when they were validated using the MCC algorithm of CYTOHUBBA, another plug-in for hub gene extraction (Fig. 3e). These results confirmed the reliability of the extracted hub genes.

Fig. 3.

The protein-protein interaction (PPI) network and its key modules extracted from the DEGs. (a) PPI network of DEGs using the STRING website. (b) Expression level of 13 hub genes in the four GEO datasets. A darker color represents a higher level of expression. (c,d) Key modules extracted by the MCODE plug-in with a degree cutoff $>$2, a node score cutoff $>$0.2, a k-score of 2, and maximum depth of 100. (e) Hub genes extracted using the MCC algorithm of CYTOHUBBA. A longer bar represents a higher score. DEG, differentially expressed gene; GEO, gene expression omnibus; PPI, protein-protein interaction.

Utilizing the information of the SIGNOR3.0 database, an open-source signaling network that manually stores annotated causal relationships between human proteins, biologically relevant chemicals, stimuli, and phenotypes, TPX2, KIF11, PTTG1, and another non-modular gene, CDK1, were extracted from 13 hub genes, forming two possible signaling pathways (Fig. 4).

Fig. 4.

Signaling pathways constructed by SIGNOR3.0. Searching for the 13 hub genes in the SIGNOR3.0 database. Two signaling pathways containing three of the hub genes were automatically generated for further analysis.

One of them comprised the CDK1, PTTG1, and mitogen-activated protein kinases/extracellular signal-regulated kinase (MAPK-ERK) signaling pathway. PTTG1 has been shown to be highly expressed in various neoplastic diseases, such as hepatocellular carcinoma, cholangiocarcinoma, and breast cancer [21, 22, 23], and is associated with tumor invasion and metastasis [24]. Both CDK1 and MAPK are kinases, whereas PTTG1 has sites for phosphorylation modification (Fig. 5a). However, activation of PTTG1 through phosphorylation by MAPK has not been confirmed in vivo, and there is a lack of studies on the PTTG1-MAPK-ERK-related pathway in cSCC.

Fig. 5.

Molecular structure and post-translational modification sites of PTTG1 (a), TPX2 (b), and KIF11 (c) according to cBioPortal.

Another pathway contained three genes, TPX2, KIF11, and CDK1, of which TPX2 and KIF11 are involved in the biological processes of spindle organization, mitotic spindle pole, and mitosis in the GO functional annotation. TPX2 includes three major functional domains: the N-terminal Aurora A kinase binding domain, the nuclear localization signal (NLS) structural domain, and the TPX2 structural domain (Fig. 5b). Among them, the TPX2 structural domain mediates the localization of Xklp2 to the mitotic spindle in Xenopus Leavis [25]. KIF11, a member of the Eg2 (KIF15) homolog family, is homologous to Xklp2. Therefore, it is reasonable to speculate that TPX2 interacts with KIF11 through the TPX2 domain and affects the function of the mitotic spindle, thereby leading to chromosomal instability [26], which is an important event in tumorigenesis. Regarding the role of CDK1 in the pathway, both TPX2 and KIF11 have a large number of sites for phosphorylation modification (Fig. 5b,c), especially TPX2, whose phosphorylation has a considerable effect on its function [27].

TPX2, KIF11, PTTG1, and CDK1 were consistently highly expressed in the four selected GEO datasets (Fig. 3b). Owing to the lack of a database involving prognosis-related clinical information for cSCC, we validated the expression levels and prognostic relevance of these four genes in a variety of squamous cell carcinomas of other tissues and cutaneous malignant melanoma that was homozygous for cSCC (Fig. 6a–h). HNSC, which refers to squamous cell carcinoma occurring mostly in the mucosal areas of the head and neck, such as the tongue, pharynx, and larynx, is broadly classified as squamous cell carcinoma of the skin. TCGA-based analyses showed that these four genes were highly expressed in most squamous cell carcinomas, as well as cutaneous melanoma (Fig. 6a–d). Moreover, elevated expression of these genes except for PTTG1 was correlated with poorer overall survival in HNSC (p 0.1; Fig. 6e–h). The HPA database was used to verify histological expression levels of TPX2 and KIF11. KIF11 was significantly more expressed in cSCC tissues than in normal tissues, whereas the difference in TPX2 staining intensity between cSCC and the adjacent normal skin tissue was not significant (Fig. 7a). We verified the expression difference of mRNA in the cell lines (Fig. 8d), whereas higher protein expression level of only TPX2 was observed in A431 cells (Fig. 8a,b).

Fig. 6.

Association of overexpressed core differentially expressed genes in different types of carcinomas with poor prognosis. Expression levels of PTTG1 (a), KIF11 (b), TPX2 (c), and CDK1 (d) in four carcinoma types (HNSC, SKCM, ESCA, and LUSC). *p 0.05. (e–h) Overall survival curves of patients with high and low expression of PTTG1 (e), KIF11 (f), TPX2 (g), and CDK1 (h) in HNSC samples. p 0.1. ESCA, esophageal carcinoma; HSNC, head and neck squamous cell carcinoma; KIFF, kinesin family member 11; LUSC, lung squamous cell carcinoma; SKCM, skin cutaneous melanoma; TPX2, targeting protein for Xklp2.

Fig. 7.

Histological expression levels and mutation likelihood of TPX2 and KIF11 according to the HPA database and cBioPortal. (a) Histological expression levels of KIF11 in the HPA. (b) Mutation rates of TPX2, KIF11, and PTTG1 in cSCC. (c,d) Mutation-prognosis curves of TPX2 and KIF11. cSCC, cutaneous squamous cell carcinoma; HPA, Human Protein Atlas; KIFF, kinesin family member 11; TPX2, targeting protein for Xklp2.

Fig. 8.

In vitro vertification of the inhibition ability to A431 cell line of NSC689534. (a,b) Expression levels of KIF11 (a) and TPX2 (b) in HaCaT cell line, A431 cell line and A431 cell line with treatment of NSC689534 at 0.25 µM for 24 h by Western Blot. ***p 0.001. (c) Inhibition rate of NSC689534 in A431 was assessed by CCK8 assay at 24 h and 36 h. IC50 (24 h) = 0.2400 µM, IC50 (36 h) = 0.2232 µM. (d) qPCR analysis to evaluate the expression of TPX2 and KIF11 in A431 and HaCaT cell lines.

We also queried the mutation status of these genes in patients with cSCC using cBioPortal. Except for PTTG1 with a mutation rate 1%, the mutation rates of TPX2 and KIF11 were 7% and 9%, respectively (Fig. 7b), and their mutations were associated with poorer overall survival (Fig. 7c,d).

Based on the CTD database, interactions between small molecular compounds and the gene TPX2 and KIF11 networks were visualized using the NetworkAnalyst website (Fig. 9a). In total, 21 common small molecule compounds were identified. After further analysis with the help of PubChem and published literature, we excluded compounds that were inherently toxic (e.g., troglitazone and mustard gas) and have been studied in cSCC (e.g., cyclosporine) [28, 29]. We selected the pyridine 2-carbaldehyde thiosemicarbazone (NSC689534) (Fig. 9b) as the target for further studies. Molecular docking analyses identified the potential of NSC689534 to interact with KIF11 (Fig. 9c). Although TPX2 could not dock directly with NSC689534, the TPX2-Aurora A complex, which affects spindle length in human somatic cells [30], was able to interact with NSC689534 through hydrophobic interactions (LEU139, VAL147, LEU263, and ASP274) and hydrogen bonds (ASN216 and ASP274) with strong binding ability (Fig. 9d and Table 3). Likewise, KIF11 bound NSC689534 strongly through four hydrophobic interactions and one $\mathrm{\Pi }$-stacking (parallel) interaction with a binding energy of -7.3 kcal/mol (Fig. 9c and Table 3). Then, we provided chemical insights into NSC689534 using SwissADME. Accordingly, the molecular weight of the compound was determined to be 362.45 g/mol. The analysis also implied a high level of gastrointestinal absorption, and the compound could not easily cross the blood-brain barrier, indicating that NSC689534 met the criteria required for a potential drug in terms of properties such as molecular weight, lipophilicity, and pharmacokinetics (Table 4). The drug-likeness indexes (Lipinski, Ghose, Veber, Egan, Muegge) were also all in compliance, and the bioavailability score was $>$0.1 (Table 4), suggesting good drug-likeness. Moreover, the synthetic accessibility score of NSC689534, a method to assess the ease to synthesize a compound, was 3.22. A score approaching one indicates that a compound has good synthetic accessibility. In a recent research, the anti-tumor capability of NSC689534 was identified in vitro and in vivo [31]. The treatment of NSC689534/Cu${}^{2+}$ at 5 µM for 24 or 48 h could induce apoptosis in HL60 and PC3 cells [31]. According to the analysis above, NSC689534 may be an effective inhibitor of TPX2 and KIF11 and show anti-tumor effect in cSCC.

Fig. 9.

NSC689534, a potential drug candidate targeting TPX2 and KIF11. (a) Small molecular compound-gene interaction network visualized on the NetworkAnalyst website. (b) Structure of NSC689534, KIF11, and TPX2. (c,d) Three-dimensional and two-dimensional structures of ligand (NSC689534)–receptor (KIF11, TPX2) interactions. KIFF, kinesin family member 11; TPX2, targeting protein for Xklp2.

A431 cells were cultured in 96-well plates and exposed to different concentrations of NSC689534 for 24 h or 48 h. After determining the cell viability with CCK8 reagent, the inhibition ability of NSC689534 to A431 cells was confirmed (Fig. 8c). It was possible to achieve very impressive inhibition rate even at a low concentration of 0.5 µM. Subsequently, the protein expression of TPX2 and KIF11 was compared among normal keratinocytes, A431 cells and cSCC cells treated with NSC689534. It was found that NSC689534 could reduce the expression of KIF11 but not TPX2 (Fig. 8a,b). Combining results of molecular docking analysis (Fig. 8c,d), we consumed that NSC689534 interacted with KIF11 and then reduced its expression in cSCC cells. Our results suggested that NSC689534 could inhibit the proliferation of A431 cells and reduce the expression of KIF11.