Identifying KRT20 as a Potential Key Gene in Lymphatic Metastasis of Head and Neck Squamous Cell Carcinoma.

Background: Head and neck squamous cell carcinoma (HNSCC) was the seventh most common cancer worldwide in 2018. Lymphatic metastasis (LM) is closely related to HNSCC prognosis and recurrence. However, the underlying mechanism of LM remains unclear. Therefore, this study aimed to identify the key genes in the LM of HNSCC.
Methods: We used The Cancer Genome Atlas (TCGA) to identify differentially expressed genes (DEGs) between LM and non-LM cases. A random forest model, the Search Tool for the Retrieval of Interacting Genes, Cytoscape, and cytoHubba were used to identify hub genes among DEGs, including KRT20 (Cytokeratins 20). We analyzed the survival of KRT20 in TCGA, and we overexpressed KRT20 in HNSCC cell lines to investigate its effects on migration and invasion. We also correlated the expression of KRT20 in HNSCC tissue microarrays with survival and clinicopathological features.
Results: We identified 243 DEGs-143 upregulated genes and 100 downregulated genes. Further analysis revealed that KRT20 is a potential key gene associated with LM and overall survival rates among patients with HNSCC. Overexpression of KRT20 increased the migration and invasion ability of HNSCC cell lines Tu686 and FD-LSC-1. Tissue microarray studies demonstrated an overexpression of KRT20 among N1+ patients (including N1-N3 patients). Survival analysis results and the clinicopathological features of HNSCC tissue microarrays were consistent with our analysis of TCGA. Thus, a high KRT20 expression level might suggest an adverse HNSCC prognosis. Our gene set enrichment analysis showed that KRT20 participates in many metabolic pathways, including those related to tumorigenesis and cancer development. Conclusions: We propose that KRT20 may be a key gene in HNSCC with LM.

KRT20 is a potential key gene associated with lymphatic metastasis and overall survival rates among HNSCC patients.

The overexpression of KRT20 increased the migration and invasion ability in HNSCC cells.

GSEA analysis revealed that KRT20 participates in many compounds’ metabolic pathways.

Introduction

Head and neck squamous cell carcinoma (HNSCC) was the seventh most common cancer worldwide in 2018 (with 890 000 new cases and 450 000 related deaths). Cervical lymphatic metastasis (LM) is among the leading causes of HNSCC recurrence and fatality. HNSCC has often been associated with lymphatic—rather than hematogenous—metastasis, and patients with LM require neck lymph node dissection and radiotherapy. Although multitherapy has significantly improved the prognosis of patients with locally advanced HNSCC, the emergence of local LM and distant metastases remains a prognostic challenge. Previous studies have shown that the gradual metastasis of primary tumor cells destroys the association between adjacent tumors and stromal cells. Metastasis also destroys blood and lymphatic vessels, leading to tumor growth after implantation or deposition into regional or distant tissues. However, the mechanism driving LM remains unclear. Therefore, an exploration of the key LM genes and the potential biomarkers of HNSCC is urgently needed.

Bioinformatics analysis has become a powerful tool for identifying potential key genes in cancer. Machine learning, such as the random forest method, can improve the classification performance of transcriptome data screening. One of this method's primary merits is its feature selection. Random forest has shown great potential to identify key therapeutic targets. Further, The Cancer Genome Atlas (TCGA) is a landmark cancer genome project that has identified the molecular features of more than 20 000 primary cancers.[5,6] The creation of TCGA and the ubiquity of high-throughput sequencing data have enhanced the sharing of genome datasets related to cancer research. In the current study, we aimed to screen LM and non-LM HNSCC data for differentially expressed genes (DEGs) and used a random forest model and a protein–protein interaction (PPI) network to identify the hub genes that affect LM among patients with HNSCC.

Materials and Methods

Sources of Expression Profiles and Clinical Data

We obtained a dataset containing the mRNA expression profiles of paired samples of HNSCC tissue. As of October 29, 2020, this dataset—which we acquired from TCGA (https://portal.gdc.cancer.gov/)—contained 298 HNSCC tissue samples, including 135 N0 samples and 163 N1+ samples (including N1-N3 patients). All of this study's sequencing procedures were performed using the Illumina HiSeq-Counts platform. All data used were obtained from TCGA. Given that the data's original deposition adhered to TCGA's ethical guidelines.

DEG Identification

We identified DEGs using the DESeq2 package. DEGs were identified using the following criteria: an absolute log2fold change (FC) > 1.0 and P < .05.

Functional and Pathway Enrichment Analyses

We performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using the DAVID 6.8 tool (https://david.ncifcrf.gov/tools.jsp). This study's cutoff for significance was set to P < .05.

PPI Network Construction and the Random Forest Model

We used the Search Tool for the Retrieval of Interacting Genes (STRING) database to analyze the PPI network of the DEGs. These results were then visually represented using Cytoscape software. The cutoff criterion for the combined score was >0.4. Subsequently, we used a cytoHubba plug-in in Cytoscape to further screen candidate genes, applying a degree algorithm with a screening criterion degree ≥ 6. We used a random forest model (ntree = 500) for DEGs with a random forest package, ranked the DEGs by their MeanDecreaseGini scores, and selected the top 30 genes. This process was repeated 3 times to verify its accuracy.

Survival Analysis of Hub Genes

We used GEPIA 2 (http://gepia2.cancer-pku.cn/#survival), a tool based on TCGA,[9,10] to analyze survival. The top 50% of the samples with higher expression levels were classified as having high expression, while the remaining 50% of the samples were classified as having low expression. We then used the Kaplan–Meier estimator with a log rank test to compare the prognostic differences between the 2 groups. The hazard ratios (HRs) are shown, along with their 95% confidence intervals and log-rank P values.

Tissue Microarray Immunofluorescence

Referring to our previous studies,[11,12] we conducted an immunofluorescence analysis of paired tissues from 68 HNSCC patients in our cohort. A tissue microarray was obtained from 80 patients (160 cores) diagnosed with HNSCC after partial laryngectomy or total laryngectomy at Eye & ENT Hospital, Fudan University. None of these patients had taken any treatment before surgery. Several cores were damaged during preparation process. Therefore, this research included 68 matched HNSCC patients. All N1+ diseases were determined by postoperative pathological examination. We used cytokeratin 20 (ab76126, abcam, 1:200 dilution) for this analysis. We also conducted a survival analysis of the hub gene in human tissue microarrays. The patients were divided into 2 groups according to their median KRT20 expression. A Kaplan–Meier estimator with a log-rank test was used to examine the prognostic difference between the 2 groups.

Statistical Analysis of Clinicopathological Features of Tissue Microarrays

We analyzed and visually represented the statistical data using GraphPad Prism 5.0 (GraphPad Software) and SPSS 22.0 software. An χ2 test was used to analyze the correlations between KRT20 expression levels and the clinicopathological features of patients with HNSCC.

Cell Culture

This study used 2 HNSCC cell lines (Tu686 and FD-LSC-1). Tu686 was provided by the China Academy of Sciences Cell Bank. FD-LSC-1 was used in our previous study. Both were laryngeal squamous cancer cell lines. Bronchial epithelial cell growth medium (BEGM, Lonza) was used to culture FD-LSC-1, at 37 °C in a humidified 5% CO2 atmosphere. Tu686 was cultured in RPMI-1640 (HyClone) supplemented with 1% (v/v) penicillin-streptomycin (Genom Biotechnology) and 10% fetal bovine serum. All cells were tested for mycoplasma contamination.

qPCR and Western Blotting

Total RNA was isolated with TRIzol reagent (Invitrogen, Thermo Fisher Scientific) and reversed-transcribed using an Evo M-MLV RT Kit (AG11728, Accurate Biology, China). The ABI 7500 Real-Time PCR System was used to perform qRT-PCR with the SYBR Green Premix Pro Taq HS qPCR Kit (AG11717, Accurate Biology, China) using primer sets synthesized by Sangon Biotech (Shanghai). The sequences of KRT20 are as follows: forward, 5′-GGTGAACTATGGGAGCGATCT-3′, and reverse, 5′-CTAGACGGTCATTTAGGTTCTGC-3′. The sequences of GAPDH are as follows: forward, 5′-TGTAGTTGAGGTCAATGAAGGG-3′, and reverse, 5′-ACATCGCTCAGACACCATG-3′.

The protein concentrations were determined using a BCA assay kit (Beyotime). Cytokeratin 20 (ab76126, abcam, 1:1000) was mixed with 5% BSA. Using standard techniques, polyvinylidene difluoride membranes were incubated with primary and secondary antibodies. Immunodetection was accomplished through enhanced chemiluminescence. Blots were digitally imaged for chemiluminescence using ECL solution (Bio-Rad) according to the manufacturer's instructions (Bio-Rad), allowing saturation to be checked. Overall, the emitted photons were quantified for each band, particularly for homogeneously loaded loading controls. The protocol was similar to that described in our previous studies.

Cell Transfection

For stable cell transfection, the KRT20 overexpression (OE) lentivirus (PGMLV-CMV-KRT20-eGFP-3Flag-PGK, Puromycin Lentivirus, Genomeditech) was mixed with serum-free RMPI and BEGM containing 1 μg/mL polybrene (Yeasen Biotechnology). After 48 h of incubation, 5 μg/mL puromycin (Sigma) was added to the culture media, and the concentration was gradually increased up to 10 μg/mL. The process lasted for 1 week to select the stably transfected cells.

Transwell Migration/Invasion Assays

We seeded 2 × 104 cells in the upper chamber of the 8 μm transwell inserts (BD Biosciences) for the transwell/migration assay with 100 μL of serum-free medium. The lower chamber was filled with a medium (50 μL) containing 10% bovine serum albumin. After 36 h of incubation, cells in the upper chamber were carefully removed. The cells that adhered to the membrane were fixed in paraformaldehyde for 15 min before being stained with 1% crystal violet (KeyGEN Biotech) for 20 min. In invasion assays, the upper chamber was coated with 50 μL of Matrigel (BD Biosciences) diluted 1:8 with serum-free medium before seeding 2 × 105 cells in 100 μL of serum-free medium. The rest of the procedure was similar to the transwell migration protocol used in our previous studies.

Gene Set Enrichment Analysis of the Whole-genome Dataset and KRT20 Expression Levels

To further examine the potential mechanism underlying diverse KRT20 expression groups in HNSCC tumor tissues, we performed a gene set enrichment analysis (GSEA). Results that met the following criteria were regarded as statistically significant: nominal P < .05, normalized enrichment score >1, and false discovery rate  < 0.3.

Results

Differentially Expressed Gene Identification and Bioinformatics Analysis

After data preprocessing and normalization using the DESeq package, we identified 243 DEGs—comprising 143 upregulated genes and 100 downregulated genes, as shown in Figure 1A and B. We conducted KEGG and GO pathway enrichment analyses using the clusterProfiler R package. To examine their cellular components, we enriched the DEGs in cornified envelopes, intermediate filaments, and extracellular space. To assess their molecular function, we also enriched the DEGs with peptide crossing-linking, keratinocyte differentiation, and keratinization. KEGG pathway clustering revealed that the DEGs were mostly enriched by serotonergic synapses, tyrosine metabolism, and salivary secretion (Figure 1C). The DEGs’ top 30 GO terms are shown in Figure 1D. Our cluster analysis of biological processes suggested that the DEGs were significantly enriched by “structural molecule activity,” “serine-type endopeptidase inhibitor activity,” and “structural constituent of epidermis.”

Figure 1.

Bioinformatic analysis of differentially expressed genes. (A) Volcano plot. (B) Heatmap. (C) Kyoto Encyclopedia of Genes and Genomes analysis. (D) Gene ontology analysis. (E) PPI network.

PPI Network Construction and Random Forest Analysis of Differentially Expressed Genes

Using STRING analysis, we input the DEGs into the PPI network, including 87 nodes and 128 sides (Figure 1E). The node indicates the protein involved in the string network, and the side indicates the direct interaction between different proteins. The details of the 87 nodes and 128 sides can be seen in Figure 1E. We then used the plugin cytoHubba to rank the nodes by their network capabilities. We also used the maximal clique centrality method to identify the top 20 genes (Table 1), including LOR, LCE3E, and KRT2, which may play important roles in HNSCC. To further identify hub genes, we used random forest packages in R to identify hub genes. The top 30 hub genes are shown in Table 2.

Table 1.

Hub Genes Identified by cytoHubba.

Rank	Name	Score
1	LOR	70
2	LCE3E	48
3	KRT2	42
4	KPRP	31
5	LCE2C	30
5	LCE2A	30
7	KRT20	26
8	LCE2D	25
9	KRT38	24
9	KRT36	24
9	KRT3	24
12	FLG2	21
12	ALB	21
14	FLG	20
15	ORM1	9
16	SLC32A1	8
16	ALOX12B	8
18	NR0B1	6
18	FGF4	6
18	FGF19	6

Table 2.

Hub Genes Identified by Random Forest Model.

GENE	MeanDecreaseGini
OTOP3	1.69517091
ACSBG1	1.62439616
OTOP2	1.5067028
GAPDHP27	1.44413024
LAMB4	1.23335332
TGM3	1.21570619
MUCL3	1.21060149
RPL7P38	1.17705076
ADH1B	1.13261751
MASP1	1.12515731
MYEOV	1.04466863
NWD2	1.04155006
DNAJB6P3	1.04024268
FRG2B	0.97080214
HTR3A	0.96627733
ROR1-AS1	0.94997745
SPINK7	0.89296173
ALOX12B	0.85848232
FRMD1	0.83678373
DMBT1	0.82331013
ALOX15B	0.81460052
LINC00452	0.80704212
PLA2G4D	0.78506928
TBC1D3P1	0.78415146
COL19A1	0.77486406
H2AC5P	0.7650745
ST8SIA6-AS1	0.75624077
CHCHD2P7	0.72569256
KRT20	0.72414612
KRT3	0.70645376

Identifying KRT20 as the Hub Gene

We used Venn plots to select intersect genes. KRT20, KRT3, and ALOX12B were identified as key genes (Figure 2A). We then conducted survival analysis of KRT20, KRT3, and ALOX12B (Figure 2B-D). For all HNSCC cases in TCGA, the results of this analysis revealed that high expression levels of KRT20 were associated with poorer overall survival rates among patients with HNSCC; this was not the case for the high expression levels of KRT3 or ALOX12B.

Figure 2.

Analysis of head and neck squamous cell carcinoma tissue microarrays. (A) Venn plot of hub genes from cytoHubba and the Random Forest Model. (B) Kaplan–Meier survival analysis of KRT3 in The Cancer Genome Atlas. (C) Kaplan–Meier survival analysis of ALOX12B in The Cancer Genome Atlas. (D) Kaplan–Meier survival analysis of KRT20 in The Cancer Genome Atlas.

KRT20 Expression is Higher Among HNSCC Patients with N1+, Indicating Poorer Overall Survival Rates

A total of 68 patients were randomly selected from our database and analyzed. KRT20 immunofluorescence staining showed higher KRT20 expression among patients with N1 + HNSCC compared to patients with N0 HNSCC (P = .0373) (Figure 3A and B). Our survival analysis showed that the overexpression of KRT20 and the N1+ stage were related to poorer overall survival rates among patients with HNSCC (P = .0203 and P = .0359, respectively) (Figure 3C and D). Further analysis revealed that KRT20 expression levels were related to cervical lymph node metastasis (Table 3).

Figure 3.

KRT20 expression is higher among N1+ head and neck squamous cell carcinoma (HNSCC) patients, indicating poorer overall survival rates. (A) Immunofluorescence analysis of KRT20 (green) and DAPI (blue) among patients with N0 HNSCC and N1 +  HNSCC tissues. (B) KRT20's relative expression level was significantly higher among patients with N1 +  HNSCC compared to patients with N0 HNSCC (P < .05). (C) Kaplan–Meier survival analysis of KRT20 among the tissue microarrays of patients with HNSCC. (D) Kaplan–Meier survival analysis of patients with N0 and N1 +  HNSCC.

Table 3.

Relationship Between KRT20 and Clinicopathological Features of HNSCC Patients.

Variable	n	KRT20 LOW (n = 33)	KRT20 HIGH (n = 35)	P-value
Age(years)				.260
≤55	8	2	6
>55	60	31	29
Gender				1.000
Female	0	0	0
Male	68	33	35
Smoking				.884
Yes	18	9	9
No	50	24	26
Drinking				.474
Yes	38	20	18
No	30	13	17
Hypertension				.059
Yes	49	20	29
No	19	13	6
Diabetes				.735
Yes	58	29	29
No	10	4	6
T-stage				.586
T1-2	18	10	8
T3-4	50	23	27
N-stage				.015
N0	34	22	12
N1+	34	11	23
M-stage				1.00
M0	68	32	36
M1	0	0	0
Location				.393
Supraglottic	22	12	10
Glottic	28	15	13
Subglottic	1	1	0
Pyriform sinus	13	4	9
Postcricoid region	2	1	1
Posterior pharyngeal	2	0	2
Lymphovascular Invasion				.590
Yes	3	1	2
No	65	32	33
Extracapsular spread				.046
Yes	11	31	26
No	57	2	9
Resection margins				.966
Positive	2	1	1
Negative	66	32	34
Perineural invasion				.258
Yes	1	1	0
No	67	32	35

Bold values represent p < 0.05.

We separated the 68 HNSCC patients into 2 groups—KRT20 high expression and KRT20 low expression—according to median KRT20 expression (the KRT20 high expression group was defined as relative fluorescence intensity >1.25). The analysis of KRT20 expression and clinical characters of HNSCC patients revealed that KRT20 expression was related to N stage (P = .015) and extracapsular spread (P = .046). In 2015, Lewis et al defined extracapsular spread as tumor cells growing outside the lymph node envelope, perilymph node fatty infiltration, or tumor cells located outside the lymph node. All patients with extracapsular spread in our study were identified by pathological examination.

Overexpression of KRT20 Promotes Migration and Invasion in HNSCC Cells

Based on a bioinformatics analysis of KRT20, we investigated the ability of KRT20 in HNSCC cell lines. KRT20 expressing lentiviruses were successfully transfected into HNSCC cell lines (Figure 4A and B). We performed transwell migration/invasion assays using Tu686 and FD-LSC-1 cells. The results of the transwell migration/invasion assays showed that more cells migrated or invaded through the chamber in the KRT20 overexpression (KRT20 OE group) group compared to the vector group (KRT20 vector group) (Figure 4C-F).

Figure 4.

KRT20 overexpression promotes migration and invasion in HNSCC cells. (A) The results of the quantitative real-time (qRT)-PCR analysis showed the mRNA levels of KRT20 expression in vector control (Vector) and overexpression (OE) groups. (B) western blotting showed the protein level of KRT20 in Vector and OE groups. (C) Representative fields in transwell migration assays of KRT20 Vector and OE groups. (D) Quantification results of the migration of cells per field. (E) Representative fields in transwell invasion assays of KRT20 Vector and OE groups. (F) Quantification results of the invasion of cells per field.

KRT20 is Related to Important Biological Processes

To further explore the molecular mechanism of KRT20 in HNSCC, we used the GSEA method to analyze a whole-genome dataset of HNSCC tumor tissues across different KRT20 expression levels. This GSEA analysis revealed that drug metabolism cytochrome P450, glycerophospholipid metabolism, glutathione metabolism, retinol metabolism, lysine degradation, valine leucine and isoleucine degradation, and steroid hormone biosynthesis were significantly associated with high KRT20 expression levels among patients with HNSCC (Figure 5).

Figure 5.

Gene set enrichment analysis was performed to further screen for the significant pathway between the different expression-level groups of KRT20 among patients with HNSCC. P < .05 was considered significant.

Discussion

In this study, we screened a total of 243 DEGs, comprising 143 upregulated genes and 100 downregulated genes, using a PPI network, and constructed a random forest model. We selected KRT20 as a potentially key gene related to LM and the prognosis of HNSCC. Our analysis of HNSCC tissue microarrays was consistent with the results of our bioinformatics analysis, both indicating that KRT20 is highly expressed in N1 + HNSCC patients and predicts a poor prognosis. Interestingly, the overall survival curve of TCGA was significantly different around 50 months, but the 2 survival curves overlapped around 100 months and seemed similar around this time point. HNSCC is a cancer of adults, with a median age at diagnosis of 66 years for HPV-negative HNSCC, 53 years for HPV-positive HNSCC, and 50 years for EBV-positive HNSCC. For elderly patients, many factors affect their overall survival, such as cardiovascular disease, and brain disease, which may lead to this phenomenon.

Cytokeratins are filament cytoskeleton proteins found in all cells. Keratins are the most prevalent intermediate filament category found in epithelial cells, and they are expressed differently depending on the cell type. Tonofilaments are keratin filaments that braid the nucleus, span the cytoplasm, and adhere to the cytoplasmic plaques of epithelial cell–cell junctions (desmosomes) in epithelia. Keratins are an integral part of the stability continuum from single cells through tissue formation, and they play a key role in the integrity and mechanical stability of both single epithelial cells and epithelial tissues. Several studies have provided evidence supporting the role of active keratin in cancer cell invasion and metastasis. The transfection of CK8 and CK18 in mouse L cells, which leads to the formation of keratin filaments and is related to higher migration and invasion ability, suggests that keratin may affect the shape and migration of cells through interaction with the extracellular environment.[16,17]

KRT20, also known as keratin 20, is a low-molecular-weight cytokeratin encoded by the KRT20 gene and located on chromosome 17q21.2. KRT20 is differentially expressed in tumors. In the context of urothelial bladder cancer, strong evidence suggests the diagnostic and prognostic value of KRT20. Moreover, molecular techniques have demonstrated higher KRT20 expression levels in bladder cancer samples compared to non-neoplastic bladder tissue. In 2010, Ye et al found that the staining intensity of KRT20 in cancer tissues was significantly correlated with bladder carcinoma grades, distant metastasis, and TNM grades. A previous study suggested similar results for colorectal tumors: KRT20 expression levels were significantly higher in colorectal tumors than in normal mucosa, and its expression was associated with colorectal cancer recurrence and median survival rates among the study's patient population. However, several studies have also indicated different roles of KRT20 in tumor progression. Researchers have found that tumors characterized by the complete absence of KRT20 expression are very poorly differentiated and contain high percentages of Ki67 + cells, contrary to our findings. According to the mentioned study, in colorectal cancer, KRT20 is expressed in enterocytes. However, our immunofluorescence results showed that KRT20 was mainly expressed in the intercellular matrix (Figure 3A). We supposed that KRT20 is involved in different processes of tumor development. More research is needed to elucidate the functions and roles of KRT20.

Several studies have also focused on biomarkers in the LM of HNSCC. In a study of a tissue cohort that included 389 primary HNSCCs, the researchers found that high ecotropic viral integration site 1 (EVI1) expression in the primary tumor correlated with lymph node metastasis, although EVI1 expression did not predict poor survival in HNSCC. Proteomic analysis of 32 HNSCC tissue samples in a 2018 study showed that several proteins may be involved in the process of LM, especially annexin A1, fatty acid-binding protein E-FABP, keratin, and type I cytoskeletal 13. Similar to our results, they found that KRT13 and KRT6A expression levels were higher in HNSCC with LM. They also revealed low expression levels of KRT1 and KRT9 in HNSCC without LM.

The role of KRT20 in HNSCC remains unclear, especially in the LM context. Through genome-wide co-expression analysis, we found that KRT20 may participate in drug metabolism-related cytochrome P450, glycerophospholipid metabolism, glutathione metabolism, retinol metabolism, lysine degradation, valine leucine and isoleucine degradation, and steroid hormone biosynthesis. The major function of cytochrome P450 is the metabolism of foreign compounds. Previous studies have demonstrated the involvement of cytochrome proteins in the development and drug resistance of various tumors. Furthermore, CYP1A1, a member of the cytochrome family, is regarded as a typical lung cancer biomarker. Glycerophospholipid metabolism plays an important role in tumor cell survival, and it affects fundamental cellular processes, including signal transduction and gene expression. Highly proliferating cancer cells require the synthesis of new fatty acids to provide a constant glycerophospholipid supply. Glutamyl transferase (GGT) is a membrane-bound enzyme that catalyzes the breakdown of extracellular GSH, allowing the creation of the constituent glutamate and cysteine required for intracellular GSH synthesis to take place. Under oxidative stress, GGT levels are considerably increased, particularly in highly metabolic cancer cells.[25,26]

The results of the current study indicated that KRT20 was highly expressed among patients with N1 + HNSCC. Moreover, KRT20 overexpression was found to be related to poorer overall survival rates among patients with HNSCC. These results were verified in our clinical sample database. We also showed that overexpression of KRT20 promoted migration and invasion in HNSCC cell lines. Accordingly, we propose that KRT20 is a potential biomarker involved in the LM of HNSCC.

However, this study has several limitations. First, it was based on TCGA database. Although the results were verified in our clinical samples containing 68 HNSCC patients, the sample size was insufficient. Importantly, this study's results must be validated using external databases. Second, our study was verified in cell experiments but not with in vivo experiments. Further, the proposed molecular mechanism and its genotype–phenotype correlation were not investigated in corresponding in vitro and in vivo experiments. Thus, the molecular mechanism of KRT20 in the LM of HNSCC remains unclear, highlighting the need for further studies.

Conclusion

Our study showed that high KRT20 expression levels among patients with HNSCC were associated with both LM and unfavorable prognoses. KRT20 may also promote the ability of migration and invasion in HNSCC cells. Through genome-wide functional enrichment analysis, we found that KRT20 may be involved in the following biological processes and pathways in HNSCC: drug metabolism-related cytochrome P450, glycerophospholipid metabolism, glutathione metabolism, retinol metabolism, lysine degradation, valine leucine and isoleucine degradation, and steroid hormone biosynthesis. However, our results must be verified by future studies.