This article is mentioned in:

Introduction

Head and neck squamous cell carcinoma is the sixth most common type of tumor worldwide, of which oral squamous cell carcinoma (OSCC) is the most frequently occurring (1,2). Every year, nearly 300,000 new cases of OSCC and ~140,000 associated deaths occur worldwide. The 5-year survival rate of patients with OSCC is ~50% (3–5). Maxillofacial defects and the loss of food and language functions that may occur in cases of OSCC severely affect the physiology and psychology of patients (6). OSCC is characterized by occult onset with high aggressiveness and invasiveness (7). Numerous patients with OSCC receive treatment in the middle or late clinical stages of the disease. Histopathological biopsies remain the primary diagnostic method for OSCC; however, there are currently no precise biomarkers for the diagnosis and prognosis of OSCC.

Biomarkers are diagnostic and prognostic tools comprising laboratory indicators associated with diagnosis or outcome of a disease (8). High-throughput sequencing technology has revealed a variety of genes that are associated with the early diagnosis, treatment and prognosis of OSCC (9). Various databases collect sequencing data and make it available to researchers; these include The Cancer Genome Atlas (TCGA) database, which has molecularly characterized >20,000 primary cancers and matched normal samples spanning 33 types of cancer (10). Differentially expressed genes (DEGs) between primary OSCCs and matched normal samples have been screened in a previous study using a dataset from TCGA (11). Furthermore, several bioinformatics tools, including Gene Ontology (GO) Resource, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Set Variation Analysis (GSVA), have been used in other studies to investigate the functions, pathways and associations with survival, DNA repair and immunocyte infiltration of DEGs (12,13).

In the present study, several notable DEGs in OSCC were identified using data from TCGA and Gene Expression Omnibus (GEO) databases. A series of bioinformatics methods and tools were used to analyze the expression of hub genes in OSCC and their association with the tumor immune microenvironment, immune checkpoints and DNA repair genes.

Materials and methods

Data collection

Gene expression profiles for human OSCC and normal oral mucosa and clinical information were obtained from TCGA database (portal.gdc.cancer.gov/). The respiratory module was used to search RNA sequencing (RNA-seq) data. The search parameters used were: Primary site: ‘Base of tongue’, ‘floor of mouth’ and ‘other and ill-defined sites in lip, oral cavity and pharynx’; Program: TCGA; Disease type: Squamous cell neoplasms; Data category: Biospecimen; Experimental strategy: RNA-seq. The criteria used when screening the data for download were as follows: Samples from oral cancer sites (alveolar ridge, base of tongue, buccal mucosa, floor of mouth, hard palate, oral cavity and oral tongue) were included; samples of non-oral cavity cancer sites (hypopharynx, larynx, lip, oropharynx and tonsil) were excluded. The mRNA expression data and clinical information of the patients, as well as chart, manifest, metadata, clinical and other files associated with the samples were downloaded using the GDC Data Transfer Tool (gdc.cancer.gov/access-data/gdc-data-transfer-tool). Finally, 361 patient samples, comprising 329 tumor samples and 32 controls with complete clinical data were included in the study.

Identification of DEGs

The DEseq2 (version 1.26.0) package in R software (version 3.6.3) (Microsoft, WA, USA) was used to identify DEGs. The original RNA-seq data was corrected by normalization to transcripts per million, using the limma package in R (version 3.6.3; bioconductor.org/packages/release/bioc/html/limma.html) for data filling, merging, correction and matrix fusion. The Ensembl gene IDs were transformed to gene names by human gene annotation using the GENCODE website (https://www.gencodegenes.org/). When identifying the DEGs, an absolute log2 fold change (|log2FC|)>2 and adjusted P<0.05 were set as the cut-off criteria. Volcano plots were generated using the R package ggplot2 (version 3.3.3). Heat maps were generated using the ComplexHeatmap R package (version 2.2.0). The expression of DKK1 in OSCC and normal tissue was analyzed based on data from TCGA database using he R package DESeq2. Two OSCC GEO Dataset (GSE3524 and GSE37991) (14,15) were also selected for confirming the expression of DKK1 in OSCC and normal tissues. The R package DESeq2 was used for differential expression analysis.

GO and KEGG enrichment analysis of DEGs

To explore the potential functions of the DEGs, functional enrichment analysis was performed. A conversion package (org.Hs.eg.db; Version 3.8) from Bioconductor was used to annotate the DEGs. The GO tool (GO;http://www.geneontology.org/) was used to categorize the genes with regard to molecular function (MF), biological pathway (BP) and cellular component (CC). GO enrichment analyses were performed using the clusterProfiler package module (version 3.14.3) in R, and the top four results of each group were identified. KEGG pathway enrichment analysis was also performed using clusterProfiler. Terms with P<0.05 were considered to be statistically significantly different and to meet the criteria and thresholds for enriched pathways.

Key gene screening

The patient population was also screened for potential prognostic genes affecting overall survival (OS) in OSCC. The patients were divided into high and low expression groups based on the median expression of the DEGs. Following Kaplan-Meier (KM) analysis of the DEGs, several significantly expressed genes in OSCC were obtained and potential key genes were searched through the VennDiagram package in R (ggplot2 version: 3.3.3). To evaluate the interactions between proteins, the STRING (http://string-db.org) online database was used. The potentially key genes were input into the module ‘Multiple proteins’, and protein interaction data with confidence >0.4 were selected as passing the threshold used to define an interaction. The protein-protein interaction (PPI) networks were then visualized and downloaded using Cytoscape software (version 3.7.2; Institute for Systems Biology).

Survival analysis

Univariate Cox analysis was used to analyze the association between Dickkopf Wnt signaling pathway inhibitor 1 (DKK1) expression and OS for various types of cancer. The KM method was used to investigate the relationship of DKK1 expression with OS and disease-specific survival (DSS). The samples were divided into DKK1 high and low expression groups based on the median expression level of DKK1. Univariate Cox survival analysis was performed using the R package survminer (version 0.4.9), and the results were visualized using the R package ggplot2 (version 3.3.3) and forest plots.

Immunocorrelation analysis

The correlations between DKK1 expression and immune cell infiltration in OSCC were analyzed using the R package GSVA (version 1.34.0). In addition, a correlation analysis of DKK1 with immune checkpoint-associated genes in OSCC was performed using the same R package.

DNA repair gene correlation analysis

The correlation of DKK1 expression in OSCC with the expression of five mismatch repair (MMR) genes, namely epithelial cell adhesion molecule (EPCAM), MutL homolog 1 (MLH1), MutS homolog 2 (MSH2), MSH6 and PMS1 homolog 2 (PMS2), was evaluated using expression profile data from TCGA. Visual analysis was performed using ggplot 2 (version 3.3.3).

Construction of lentiviral particles

Lentiviral particles carrying a short hairpin RNA (shRNA) targeting DKK1 or a shRNA control (Lv-shDKK1 and Lv-shCon, respectively) were designed and synthesized by Shanghai Genechem. Briefly, shCon and shDKK1 were inserted into a GV248 lentiviral vector (Shanghai GeneChem). The vector (1.5 µg) carrying shCon or shDKK1 and Lenti-Easy Packaging Mix (Shanghai GeneChem) (1.5 µg) were transfected into 293T cells (#bio-12947; Biobw) with Lipo6000™ reagent (#C0526; Beyotime Institute of Biotechnology) at 37°C for 24 h, and the lentiviral particles were collected by centrifugation at 80,000 × g and 4°C for 2 h using an ultra-centrifuge (Beckman Coulter, Inc.). The shRNA sequences were as follows: sh-DKK1, 5′-CACCGCTCTCATGGACTAGAAATATTTGATATCCGATATTTCTAGTCCATGAGAGC-3′; and shCon, 5′-CACCTCTGTCAATTAGGACAAGCTTATGATATCCGTAAGCTTGTCCTAATTGACAGA-3′. Subsequent experiments were performed 24 h later.

Cell culture

The OSCC cell lines CAL-27 and SCC15 were purchased from the American Type Culture Collection. The cells were cultured in Dulbecco's modified Eagle's medium (DMEM; cat. no. BL304A; Biosharp Life Sciences) containing 10% fetal bovine serum (FBS; cat. no. AB-FBS0500; Animal Blood Ware) at 37°C with 5% CO2.

Lentiviral infection

CAL-27 and SCC15 cells were seeded into a 12-well plate at 20,000 cells per well. When a confluence of 40–50% was reached, Lv-shCon or Lv-shDKK1 at a multiplicity of infection of 10 was added to each well. After infection for 24 h, the medium containing the lentivirus was replaced with DMEM without lentivirus. At 72 h after infection, the cells infected with the lentivirus were selected using 1 µg/ml puromycin (Beyotime Institute of Biotechnology). Reverse transcription-quantitative PCR (RT-qPCR) was used to confirm the knockdown of DKK1.

RNA extraction and RT-qPCR

Cells were lysed using TRIzol® reagent (cat. no. 15596; Thermo Fisher Scientific, Inc.), and RNA was extracted from the lysates using chloroform and separated using isopropanol. The RNA concentration was determined using an ultramicro spectrophotometer (Implen GmbH). The RNA was reverse-transcribed into cDNA using a BeyoRT™ M-MuLV Reverse Transcriptase kit (cat. no. 7268; Beyotime Institute of Biotechnology) according to the manufacturer's instructions. The cDNA was then subjected to qPCR using a BeyoFast™ SYBR Green qPCR Mix kit (cat. no. 7260; Beyotime Institute of Biotechnology) according to the manufacturer's instructions following the thermocycling conditions: 95°C 10 min; (95°C, 15 sec; 60°C, 1 min) ×40 cycles. The primers used were as follows: DKK1 forward, 5′-GAATAAGTACCAGACCATTGAC-3′ and reverse, 5′-CCATTTTTGCAGTAATTCCC-3′; β-actin forward, 5′-CGGGAAATCGTGCGTGAC-3′ and reverse, 5′-CAGGCAGCTCGTAGCTCTT-3′. The expression of target genes was quantified using β-actin as the reference gene via the 2−ΔΔCq method (16).

Wound healing assay

CAL-27 and SCC-9 cells were seeded in 6-well plates at a density of 500,000 cells/well. When the cells reached 90% confluence, a wound was created by drawing a straight line through the cells using a pipette tip. The cells were then cultured in DMEM without FBS, and the wounds were observed and evaluated under an inverted light microscope (IX53, Olympus, Tokyo, Japan) at 0, 24 and 48 h. The wound healing rate was calculated as follows: (Original wound area-non-healing wound area)/original wound area ×100.

Transwell assay

The invasion capacity of the cells was evaluated using a Transwell assay. Briefly, the upper Transwell membrane was pre-coated with Matrigel (E6909; Sigma-Aldrich; Merck KGaA) at 37°C for 4 h. CAL-27 and SCC-9 cells were starved for 24 h in serum-free DMEM to stop cell proliferation. After that, 100,000 cells in 200 µl FBS-free DMEM were seeded in the upper chamber of the Transwell system, while 500 µl medium containing 10% FBS was placed in the lower chamber. After culture for 12 h at 37°C, cells crossing the membrane were stained with crystal violet (cat. no. C0121; Beyotime Institute of Biotechnology) for 3 min at 20°C. The cells were then observed and counted under a light microscope (IX53, Olympus, Tokyo, Japan).

CCK-8 assay

A CCK-8 assay was used to evaluate cell proliferation ability. Briefly, the CAL-27 and SCC-9 cells were seeded in a 96-well plate at 5,000 cells/well. After culture for 1, 24 and 48 h, cells were incubated in 10 µl CCK-8 reagents (cat. no. C0037; Beyotime Institute of Biotechnology) for 1 h at 37°C. Cell proliferation was determined by calculating the optical density at 450 nm.

Colony formation assay

The ability of the cells to form colonies was determined using 6-well plates seeded with 1,000 cells/well. The CAL-27 and SCC-9cells were incubated for 12 days, then fixed with 4% paraformaldehyde for 5 min at 20°C, washed with phosphate-buffered saline and stained with crystal violet for 5 min at 20°C. Clonogenicity was observed and the colonies were photographed using an inverted microscope. The colonies larger than 0.3 mm were calculated using imageJ software (version 1.8.3) (NIH).

Western blot analysis

Proteins were extracted from SCC-9 and CAL-27 cells by lysis using RIPA buffer (cat. no. C1053; Applygen Technologies, Inc.). Protein concentration was determined using a BCA kit (cat. no. P0010; Beyotime Institute of Biotechnology) and total proteins (30 µg/lane) were separated using 10% SDS-PAGE and then transferred to a PVDF membrane (cat. no. FFP22; Beyotime Institute of Biotechnology). The membrane was blocked with 5% fat-free milk (cat. no. P0216; Beyotime Institute of Biotechnology) at 25°C for 1 h, incubated with primary antibodies [1:1,000 diluted in Tris-buffered saline with 1% Tween 20 (TBST)] for 12 h at 4°C and then incubated with secondary antibodies (1:10,000 diluted in TBST) for 1 h at 20°C. Finally, the membranes were incubated in ECL reagent (cat. no. P0018; Beyotime Institute of Biotechnology) for 1 min at 20°C and developed using a chemiluminescence imager (Tanon Science & Technology Co., Ltd.). β-actin antibody (cat. no. sc-8432) was purchased from Santa Cruz Biotechnology, Inc. Wnt-3a antibody (cat. no. 2391), DKK1 antibody (cat. no. 4687) and β-catenin (cat. no. 8480) antibody were purchased from Cell Signaling Technology, Inc. Goat anti-rabbit HRP-secondary antibody (cat. no. A0208) and goat anti-mouse HRP-secondary antibody (cat. no. A0216) were purchased from Beyotime Institute of Biotechnology.

ELISA

Secreted DKK1 in the culture supernatant was detected using a Human DKK1 ELISA Kit (cat. no. EK0867; Wuhan Boster Biological Technology, Ltd.) following the manufacturer's instructions. Briefly, 100 µl supernatant or standards were added to the plate pre-coated with a DKK1 antibody. Subsequently, 2 µl DKK1 detector antibody was added and the plate was incubated in a second antibody with horseradish peroxidase for 1 h at 37°C. After washing the plate three times with wash solution, ABC working reagent was added. The color reaction was performed by adding 90 µl tetramethylbenzidine reagent for 20 min in the dark and terminated by adding 100 µl stopping reagent. The optical density at 450 nm was measured using a microplate reader (Molecular Devices, LLC). A standard curve was established based on the optical density of the standards and used to calculate the DKK1 concentration in the supernatant.

Cell apoptosis

Cell apoptosis was detected using an AnnexinV/PI staining kit (cat. no. APOAF-20TST; Sigma-Aldrich; Merck KGaA). Briefly, 5×105 SCC-9 or CAL-27 cells were harvested and suspended in 500 µl binding buffer. Then 5 µl Annexin V/PI staining reagent was added and mixed with the cell suspension. Following incubation for 10 min, cell apoptosis was determined using a flow cytometer (BD LSRFortessa™; BD Biosciences) and analyze using the FlowJo software (version 10.4.0) (FlowJo, OR, USA).

Statistical analysis

The difference in gene expression between the two groups was analyzed using the Wilcoxon rank-sum test, except for GSE37991 for which Welch's t-test was used. Cox regression analysis was used to determine the survival probabilities of the two groups. KM plots were drawn and the relationships between various factors and survival outcomes were analyzed using Cox regression. The correlation of DKK1 expression with immune factors and DNA repair genes in OSCC was analyzed using Spearman's correlation method. Data from the RT-qPCR, wound healing, Transwell, CCK-8 and colony formation assays are presented as the mean ± SD and were analyzed using unpaired Student's t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

DEGs in OSCC

The clinical features of 361 samples, including 329 samples of OSCC tissue and 32 samples of normal control tissue, were obtained from TCGA data portal. The total number of Gene IDs was 56,494, and the number of those satisfying the threshold |log2(FC)|>2 and P<0.01 was 3,226; the number with higher expression in the tumor group than the normal group (positive logFC) was 1,524, and the number with lower expression in the tumor group than the normal group (negative logFC) was 1,702 as determined using the DEseq2 package with |log2(FC)|>2, P<0.05 (Fig. 1A).

Figure 1. DEGs and functional enrichment analysis of OSCC. (A) Volcano plot of the DEGs. Red dots indicate significantly upregulated DEGs and blue dots indicate significantly downregulated DEGs in OSCC (absolute log2 fold change >2 and P.adj <0.01). (B) Heatmap plot of the DEGs. (C) Gene Ontology and (D) Kyoto Encyclopedia of Genes and Genomes enrichment analysis. DEGs, differentially expressed genes; OSCC, oral squamous cell carcinoma; P.adj, adjusted P-value.

The heat map in Fig. 1B shows the top 100 DEGs. These 100 DEGs were subjected to GO enrichment analysis. The DEGs were significantly enriched in the BP terms ‘muscle contraction’, ‘muscle system process’, ‘myofibril assembly’ and ‘muscle filament sliding’. In addition, the enriched MF terms comprised ‘structural constituent of muscle’, ‘extracellular matrix structural constituent’, ‘receptor ligand activity’ and ‘endopeptidase inhibitor activity’ and the enriched CC terms comprised ‘contractile fiber’, ‘myofibril’, ‘contractile fiber part’ and ‘sarcomere’ (Fig. 1C). KEGG pathway analysis revealed that the DEGs were mainly enriched in the pathways ‘salivary secretion’, ‘protein digestion and absorption’, ‘dilated cardiomyopathy’, ‘hypertrophic cardiomyopathy’, ‘cardiac muscle contraction’, ‘calcium signaling pathway’ and ‘neuroactive ligand-receptor interaction’ (Fig. 1D).

Screening the prognostic genes of OSCC

Through the univariate Cox analysis of genes in OSCC, 2,130 protein-coding genes were found to be associated with OS. Comparison of these 2,130 genes with the 3,226 DEGs identified that 135 key DEGs in OSCC were common to both datasets (Fig. 2A). The STRING online database was used to determine interactions among the proteins encoded by the key DEGs, and a PPI network was constructed (Fig. 2B). Based the PPI network, the top 20 genes according to statistical significance (10 positive and 10 negative associations) were selected, among which DKK1 showed the highest association with OSCC OS (hazard ratio=1.903; Fig. 2C).

Figure 2. Identification of DEGs associated with OS in OSCC. (A) Venn diagram of DEGs and OS-related genes in OSCC. (B) Protein-protein interaction network of key DEGs constructed using Cytoscape. (C) Forest plot of the top 20 genes associated with OS. DEGs, differentially expressed genes; OS, overall survival; OSCC, oral squamous cell carcinoma; HR, hazard ratio.

Selection and validation of hub gene signatures

The expression of DKK1 in OSCC and normal tissues was analyzed based on data from TCGA database, which showed that DKK1 expression levels were higher in OSCC tissue than in normal control tissues (Fig. 3A). To verify the expression of DKK1 in OSCC, DKK1 expression was analyzed in two GEO gene expression datasets: GSE3524 and GSE37991. In GSE3524, the expression of DKK1 in OSCC was higher than that in normal oral mucosa, and the median difference between the two groups was 1.285 (0.716–2.102), which was statistically significant (Fig. 3B); the analysis of GSE37991 provided a similar result, as the median difference between the two groups was 1.224 (0.978–1.47), which was also statistically significant (Fig. 3C). The expression of DKK1 in tumor and normal tissues was also compared in 33 tumor types based on TCGA data, as shown in Fig. 3D. The results revealed that DKK1 expression was significantly changed in 12 types of tumor, among which DKK1 was upregulated in colon adenocarcinoma, cholangiocarcinoma (CHOL), esophageal carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), liver hepatocellular carcinoma, lung squamous cell carcinoma (LUSC), stomach adenocarcinoma (STAD) and thyroid carcinoma compared with the corresponding normal tissues, and downregulated in bladder urothelial carcinoma (BLCA), kidney chromophobe, kidney renal papillary cell carcinoma and prostate adenocarcinoma (PRAD) compared with the corresponding normal tissue (Fig. 3D).

Figure 3. Expression of DKK1 in various tumor types. Expression levels of DKK1 in (A) The Cancer Genome Atlas, (B) GSE3524 and (C) GSE37991 oral squamous cell carcinoma datasets. (D) Expression level of DKK1 in 33 cancer types. *P<0.05, **P<0.01, ***P<0.001. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; TPM, transcripts per million; ns, not significant; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, lymphoid neoplasm diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma.

Prognostic analysis of DKK1 expression pan-cancer

The association of DKK1 expression with OS and DSS in the 12 types of tumor was calculated using univariate survival analysis. The results shown in Fig. 4A indicate that DKK1 was significantly associated with OS in ESCA, HNSC and STAD. KM analysis shows that low DKK1 expression was associated with poor OS prognosis in patients with ESCA (Fig. 4B), while high DKK1 expression was associated with poor OS prognosis in patients with HNSC (Fig. 4C) and STAD (Fig. 4D). DKK1 expression was also found to be significantly associated with DSS in HNSC and STAD (Fig. 4E), and KM analysis suggested that high DKK1 expression was associated with poor prognosis in patients with HNSC (Fig. 4F) and STAD (Fig. 4G).

Figure 4. Univariate survival analysis of the association between DKK1 expression and survival time in 12 tumor types. (A) Forest plot showing the relationship between DKK1 expression and OS. KM curves of high and low DKK1 expression in (B) ESCA, (C) HNSC and (D) STAD reveal a significant association with OS. (E) Forest plot showing the relationship between DKK1 expression and DSS. KM curves of high and low DKK1 expression in (F) HNSC and (G) STAD show a significant association with DSS. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; OS, overall survival; DSS, disease-specific survival; KM, Kaplan-Meier; BLCA, bladder urothelial carcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; ESCA, esophageal carcinoma; HNSC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRP, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; LUSC, lung squamous cell carcinoma; PRAD, prostate adenocarcinoma; STAD, stomach adenocarcinoma; THCA, thyroid carcinoma; HR, hazard ratio.

Correlation analysis of DKK1 and immune cells

The correlations between DKK1 expression and the infiltration of 24 immune cell types in OSCC were analyzed (Fig. 5A). Among these immune cells, T, regulatory T (TReg), plasmacytoid dendritic cell (pDC), T helper 17 (Th17), T follicular helper (TFH), cytotoxic and B cells showed a significant negative correlation with DKK1 expression. However, Th2 cells and gdT (Tgd) cells showed a significant positive correlation with DKK1 expression in OSCC (Fig. 5B-J).

Figure 5. Correlation analysis of DKK1 and immune cell infiltration. (A) Correlation of DKK1 with the infiltration of a panel of immune cells. Negative correlation of DKK1 expression with the level of infiltration of (B) T cells, (C) Treg cells, (D) pDC cells, (E) Th17 cells, (F) TFH cells, (G) cytotoxic cells and (H) B cells. Positive correlation of DKK1 expression with the level of infiltration of (I) Th2 cells and (J) Tgd cells. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; DC, dendritic cell; aDC, activated DC; iDC, inflammatory DC; pDC, plasmacytoid DC; NK, natural killer; Tcm, T central memory; Tem, T effector memory; TFH, T follicular helper; Tgd, γδT; Th, T helper; TReg, regulatory T; TPM, transcripts per million.

Correlation between DKK1 and immune checkpoints

Correlation analysis between DKK1 and genes associated with immune monitoring checkpoints in OSCC was performed using the R package GSVA. In total, 15 immune checkpoint genes were identified that correlated with DKK1 in the OSCC samples (Fig. 6). Among these genes, CD40, CD44, v-set domain containing T cell activation inhibitor 1 (VTCN1), neuropilin 1 (NRP1), programmed cell death 1 ligand 2 (PDCD1LG2), CD276, CD80 and CD86 were positively correlated with DKK1 expression, whereas TNF receptor superfamily member 18 (TNFRSR18), CD27, T cell immunoreceptor with Ig and ITIM domains (TIGIT), indoleamine 2,3-dioxygenase 2 (IDO2), CD48, CD244 and CD40 ligand (CD40LG) were negatively correlated with DKK1.

Figure 6. Correlation between DKK1 and genes associated with immune checkpoints in oral squamous cell carcinoma. DKK1 expression was positively correlated with the expression levels of (A) CD40, (B) CD44, (C) VTCN1, (D) NRP1, (E) PDCD1LG2, (F) CD276, (G) CD80 and (H) CD86 and negatively correlated with expression levels of (I) TNFRSR18, (J) CD27, (K) TIGIT, (L) IDO2, (M) CD48, (N) CD244 and (O) CD40LG. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; VTCN1, v-set domain containing T cell activation inhibitor 1; NRP1, neuropilin 1; PDCD1LG2, programmed cell death 1 ligand 2; TNFRSR18, TNF receptor superfamily member 18; TIGIT, T cell immunoreceptor with Ig and ITIM domains; IDO2, indoleamine 2,3-dioxygenase 2; CD40LG, CD40 ligand; TPM, transcripts per million.

DNA repair gene correlation analysis

The correlation of DKK1 expression in OSCC with the DNA MMR genes EPCAM, MLH1, MSH2, MSH6 and PMS2 was evaluated using expression profile data from TCGA. As shown in Fig. 7, DKK1 expression was significantly correlated with the DNA repair genes MSH2, MSH6, PMS2 and EPCAM but not with MLH1.

Figure 7. Correlation of DKK1 expression with that of five DNA mismatch repair genes in oral squamous cell carcinoma. *P<0.05; **P<0.01; ***P<0.001. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; MLH1, MutL homolog 1; MSH2/6, MutS homolog 2/6; PMS2, PMS1 homolog 2; EPCAM, epithelial cell adhesion molecule.

Knockdown of DKK1 inhibits cell proliferation, clonogenicity, migration and invasion in OSCC cells

To further elucidate the role of DKK1 in OSCC progression, its influence on cell proliferation, migration and invasion was investigated in vitro. DKK1 was knocked down in the OSCC cell lines CAL-27 and SCC-9 using lentiviral infection. As shown in Fig. 8A and B, the green fluorescence signal indicated that the SCC-9 and CAL-27 cells were infected with Lv-shCon or Lv-shDKK1. The knockdown efficiency was evaluated by RT-qPCR (Fig. 8C); the knockdown efficiencies of DKK1 in the SCC-9 and CAL-27 cells were 71.15 and 66.59%, respectively. Considering that DKK1 is a secreted protein, the secreted form of DKK1 was detected by ELISA, and was observed to decrease significantly after DKK1 knockdown (Fig. 8D). The role of DKK1 in cell proliferation, colony formation, migration and invasion was determined. The proliferation of SCC-9 and CAL-27 cells was inhibited by DKK knockdown (Fig. 8E and F), and the colony forming capacity of the cells was also reduced (Fig. 8G). Furthermore, the wound healing assay showed that the knockdown of DKK1 inhibited cell migration (Fig. 8H-K), and the Transwell assay indicated that the knockdown of DKK1 inhibited cell invasion (Fig. 8L-N). Furthermore, the expression of proteins associated with Wnt signaling downstream of DKK1 was investigated, and the results indicated that Wnt-3a and b-catenin levels increased after DKK1 knockdown (Fig. 8O).

Figure 8. Knockdown of DKK1 inhibits cell proliferation, colony formation, migration and invasion in oral squamous cell carcinoma cells. Bright-field and GFP fluorescence images of (A) SCC-9 cells and (B) CAL-27 cells with Lv-shCon and Lv-shDKK1 infection. Scale Bar=50 µm. (C) DKK1 mRNA expression and (D) DKK1 secretion by SCC-9 cells and CAL-27 cells infected with Lv-shCon and Lv-shDKK1. Cell proliferation rate of (E) SCC-9 and (F) CAL-27 cells. (G) Colony formation of SCC-9 cells and CAL-27 cells. (H) Representative wound healing images of SCC-9 cells and (I) quantification of the wound healing assay. (J) Representative wound healing images of CAL-27 cells and (K) quantification of the wound healing assay, magnification: 40×. Representative images of the Transwell invasion assay for (L) SCC-9 and (M) CAL-27 cells, magnification: 100×. (N) Cell counts of invaded SCC-9 and CAL-27 cells. (O) Western blot analysis of Wnt-3a, β-catenin and DKK1 expression in SCC-9 and CAL-27 cells infected with Lv-shCon and Lv-shDKK1. *P<0.05 **P<0.01 and ***P<0.001 for shDKK1 vs. shCon. DKK1, Dickkopf Wnt signaling pathway inhibitor 1; GFP, green fluorescent protein; Lv, lentivirus; shCon, shRNA control; shDKK1, shRNA targeting DKK1; shRNA, short hairpin RNA.

Discussion

Accurate cancer biomarkers can be used to predict the prognostic risk of patients and formulate individual therapeutic strategies (17). The present study aimed to identify the hub DEGs associated with OSCC prognosis and explore their potential contributions to biological processes and functions. Bioinformatic tools were used as they can effectively compensate for the shortcomings of sequencing data analysis and integrate various existing sequencing data to further investigate their clinical significance (18).

In the present study, data on OSCC from TCGA database were filtered and 3,226 DEGs were identified. In addition, 2,130 genes associated with OSCC OS-were identified using univariate Cox analysis. Furthermore, by examining the intersection of the DEGs with the OS-related genes, 135 key genes were identified, among which DKK1 showed the greatest prognostic association with OS.

According to the dataset from TCGA and the two GEO datasets, GSE3524 and GSE37991, DKK1 is upregulated in OSCC tissues. However, DKK1 is not upregulated in all tumor types; as TCGA database suggests, DKK1 is upregulated in several tumors, including CHOL, ESCA, LUSC and STAD, whereas it is downregulated in other tumors, including BLCA and PRAD. These results indicate that the expression of DKK1 may be tumor type-specific. Subsequently, the prognostic value of DKK1 was analyzed in various types of cancer, and the results showed that the expression of DKK1 was significantly increased in HNSC and STAD compared with normal tissue, and that the OS and DSS of the high-expression DKK1 group were significantly lower than those of the low-expression DKK1 group in HNSC and STAD, indicating the prognostic value of DKK1 in these tumors. In the in vitro experiments performed in the current study, DKK1 knockdown decreased OSCC cell migration and invasion. These results are inconsistent with those in a previous study (19), which showed that DKK1-positive cases were significantly associated with a low risk of regional lymph node metastasis, and cellular migration and invasion were negatively regulated by DKK1 knockdown. However, a report by Wang et al (20) supports the present study, as it showed that the proliferation and migration of OSCC cells were inhibited by the inhibition of DKK1. Therefore, the role of DKK1 in OSCC cell migration and invasion is unclear. In the study by Ogoshi et al, it was shown that DKK1 regulates the phosphorylation of β-catenin in cell nuclei; therefore, the authors speculated that DKK1 regulates cellular migration and invasion through the regulation of β-catenin phosphorylation in nuclei. In the present study, it was found that the knockdown of DKK1 increased the expression of β-catenin. Notably, Sa3 and H1 cell lines were used in the study of Ogoshi et al, while the SCC-4 cell line was used in the study of Wang et al and CAL-27 and SCC-9 cells were used in the present study. Different cell lines are derived from different patients, who have individual differences. Disparities in the reaction of β-catenin to DKK1 in different cells may be the reason for the conflicting roles observed for DKK1 in OSCC cell migration and invasion. Although it has been widely reported that DKK1 is associated with cancer cell migration and invasion, the function of DKK1 in different tumor types may vary. For example, the reduction of DKK1 expression in ovarian cancer cells has been shown to induce cell migration and invasion via activation of the serine/threonine-protein kinase 3/FOXO3 pathway (21). This may be due to DKK1 interacting with Wnt receptors, thereby inhibiting the classical Wnt signaling pathway and promoting tumor invasion and migration (22,23). However, DKK1 has also been shown to inhibit the invasion and migration of breast cancer cells by suppressing the β-catenin/matrix metalloproteinase 7 pathway (24). The reason that DKK1 functions differently in different tumor types requires further analysis.

The results of the CCK-8 assay in the present study indicated that the proliferation of the DKK1 knockdown cells was decreased compared with that of the control cells; therefore, it was important to avoid the effect of growth suppression on cell migration and invasion during the Transwell assay. Notably, in the CCK-8 assay, the cells were cultured in DMEM with 10% FBS. However, in the Transwell assay, the cells were cultured in serum-free DMEM for 24 h to stop cell proliferation. Therefore, during the Transwell assay, there was no marked cell proliferation. Also, duration of the Transwell assay was limited to 12 h. In addition, during cell culture, DKK1 knockdown was not observed to induce any obvious cell death when viewed under a microscope (Fig. S1A). Cell death was also analyzed using flow cytometry, and the knockdown of DKK1 was found to have no effect on cell death (Fig. S1B). Therefore, in the present study, the possibility that the inhibitory effects of shDKK1 on invasion and migration were due to cell death may be discounted. In addition, Transwell assays have been widely performed in numerous high-quality studies (25–27), without consideration of the effects of cell proliferation and cell death.

Immune cells in the tumor microenvironment play important roles in the occurrence and development of the tumor and are significantly associated with prognosis (28,29). They participate in remodeling the microenvironment, and regulating tumor progression; the tumor microenvironment affects immune cell infiltration (30), and targeting genes involved in this process is a promising strategy for tumor therapy. In a mouse model of ovarian cancer, the overexpression of DKK1 was found to decrease the infiltration of CD45+ leukocytes into the peritoneum and omentum, and to reduce the numbers of natural killer and CD8 T cells and the expression of interferon-g on activated CD8 T cells (31). In other studies, DKK1 was demonstrated to be associated with antitumor immunity and serve as a potential predictive marker and target for immunotherapy in several types of tumors, including lung adenocarcinoma and endometrial carcinoma (32–35). However, whether DKK1 has potential as a biomarker and target for OSCC is not clear. Therefore, the correlation between DKK1 expression and 24 immune cell types in OSCC was analyzed in the present study. The results showed that DKK1 was negatively correlated with the infiltration of T, TReg, pDC, Th17, TFH, cytotoxic and B cells, whereas it was positively correlated with the infiltration of Th2 and Tgd cells in OSCC. In addition, DKK1 positively correlated with the expression of several immune checkpoint genes, including CD40, CD44, VTCN1, NRP1, PDCD1LG2, CD276, CD80 and CD86, and negatively correlated with TNFRSR18, CD27, TIGIT, IDO2, CD48, CD244 and CD40LG. These findings suggest that DKK1 may affect the immune status of the tumor in OSCC via the regulation of specific immune checkpoint genes. Notably, DKK1 is an inhibitor of Wnt signaling, which can affect the tumor immune status in various types of tumor (36–38). However, studies on the relationship between Wnt signaling and the DKK1-related immune checkpoint genes screened in the present study are lacking. It may be hypothesized that the Wnt signaling pathway is responsible for the association of DKK1 with tumor immunity.

The DNA MMR system recognizes and corrects occasional DNA base mismatches in non-homologous chromosomes during DNA replication, and this correction ensures the stability and integrity of the genome (39). Deficiencies in the MMR system can lead to genetic mutations and induce tumorigenesis (40). Various genes participate in DNA MMR, including members of the MSH family, BRAF, PMS2 and EPCAM (41,42). In the present study, DKK1 expression was found to significantly correlate with the expression of the DNA repair genes MSH2, MSH6, PMS2 and EPCAM, suggesting that DKK1 participates in the regulation of DNA repair in OSCC. MSH2, MSH6, PMS2 and EPCAM have been reported to be associated with the Wnt signaling pathway (43–45). Therefore, it is important to verify whether DKK1 participates in DNA MMR through the Wnt signaling pathway.

Although a search of the literature found that several of the DNA repair genes and immune checkpoints screened in the present study are associated with the Wnt pathway (Wnt/β-catenin signaling induces DNA damage repair in ameliorating radio-resistance (46), no mechanistic studies were found. In future studies, detection of the mechanism by which DKK1/Wnt regulates the transcription of these genes will be investigated.

One shortcoming of the present study is that it predominantly involves bioinformatics analysis. Further in vivo and in vitro studies are required to verify the role of DKK1 in OSCC progression and to reveal the mechanism by which DKK1 participates in OSCC immunity, DNA MMR, cell proliferation, cell migration and invasion. Another shortcoming is that only DKK1 knockdown was performed; the DKK1 overexpression experiments that could provide a more convincing conclusion were omitted. In addition, the expression levels of Wnt3a and β-catenin were detected in cells with DKK1 knockdown. β-catenin can transfer into the nucleus to affect the transcription of target genes (47). Therefore, in addition to showing that DKK1 knockdown upregulates β-catenin expression, it would also be interesting to investigate whether β-catenin is also activated by DKK1 knockdown.

In summary, the present study explored the relationship between DKK1 and OSCC prognosis and the possible underlying mechanisms. The results demonstrated that DKK1 is upregulated in OSCC, and associated with survival, tumor immunity, DNA MMR and cell proliferation, migration and invasion. These findings indicate that DKK1 is a candidate gene for OSCC therapy and prognosis.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The study was supported by the National Natural Science Foundation of China (grant no. 81602374), the Natural Science Foundation of Shandong Province (grant no. ZR2021MH176), the China Postdoctoral Science Foundation (grant no. 2021M701538) and the Natural Science Foundation of Liaocheng People's Hospital (grant no. LYQN201903).

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the Figshare repository (https://figshare.com/search?q=10.6084%2Fm9.figshare.21671213).

Authors' contributions

ZM and GXZ designed the study. YJL and CCW performed the cell culture and bioinformatics analysis. SW performed lentiviral infection and drafted the manuscript. SXD performed the RT-qPCR experiments. YGL and DPZ performed CCK-8 and transwell assay and confirm the authenticity of the raw data. YNL performed the ELISA and western blot analysis. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References