This article is mentioned in:

Introduction

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common malignancy worldwide (1). There are ~650,000 new cases of HNSCC diagnosed every year, with >350,000 associated deaths occurring in a year (2). Several risk factors involved in the development of HNSCC have been identified, including alcohol, tobacco use and human papillomavirus (HPV) infection, the latter of which is a major independent risk factor for HNSCC (3). Compared with HPV-positive HNSCC, HPV-negative HNSCC is more aggressive and is associated with a higher rate of drug resistance (4). Similarly, patients with HPV-negative HNSCC demonstrate lower overall survival rates than those with HPV-positive HNSCC (5). The identification of novel and specific biomarkers for HPV-negative HNSCC may therefore help to improve the understanding of the mechanism that underlies the progression of HPV-negative HNSCC, as well as to improve the diagnosis and treatment of the disease.

With the progression of high-throughput sequencing technology and bioinformatics methods, additional oncogenes involved in the development of HNSCC have been identified (6,7). Weighted gene co-expression network analysis (WGCNA) is a bioinformatics method that involves calculating the associations between gene co-expression modules and clinical traits using a weighted soft threshold (8). Recently, with the use of WGCNA, a series of biomarkers for HNSCC have been identified. For example, Song et al (9) identified 16 genes involved in the immune and inflammatory response that promoted the development of HNSCC. Similarly, using WGCNA, Zhang et al (10) identified 12 genes associated with perineural invasion in patients with HNSCC. However, the molecular mechanism of each HNSCC subtype is different, and the feasibility and specificity of the majority of biomarkers obtained from previous studies are limited.

The current study combined WGCNA, differential expression gene (DEG) analysis and experimental verification to identify specific biomarkers for HPV-negative HNSCC.

Materials and methods

Clinical specimens

A total of 62 pairs of primary HNSCC tissues and adjacent non-tumor tissues (1-2 cm distance from the tumor tissues) were collected from patients at The Affiliated Hospital of Guizhou Medical University (Guiyang, China). None of the enrolled patients underwent treatment (such as chemotherapy or radiation) prior to tissue collection. HPV infection status had been determined prior to study commencement by performing in situ hybridization (ISH), as previously described (11), and P16 staining (Fig. S1). HPV-positive HNSCC was diagnosed based on ISH or P16 (also named p16INK4a) positive results, while HPV-negative HNSCC was defined based on negative ISH and P16 results. A total of 35 patients with HPV-positive HNSCC (26 tissues were ISH and P16 positive, 9 tissues were ISH positive only and none of the tissues were P16 positive only) and 27 patients with HPV-negative HNSCC were enrolled in the present study. The basic clinical characteristics (American Joint Committee on Cancer) (12) of the patients with HPV-negative and HPV-positive HNSCC are presented in Table I. The current study was approved by the Human Trait Ethics Committee of Guizhou Medical University and performed in accordance with the Declaration of Helsinki. All patients who provided samples signed written informed consent.

Table I. Basic clinical characteristics of patients with HPV-negative (n=27) and HPV-positive (n=35) HNSCC.

Table I.

Basic clinical characteristics of patients with HPV-negative (n=27) and HPV-positive (n=35) HNSCC.

Clinicopathological feature	HPV-negative HNSCC	HPV-positive HNSCC	P-value
Mean age ± SD, years	71.4±6.7	72.4±7.3	>0.05
Sex, n (%)			>0.05
Male	12 (44.4)	17 (48.5)
Female	15 (55.6)	18 (51.5)
Tumor grade (Broders), n (%)			>0.05
Grade I	12 (44.4)	16 (45.7)
Grade II	13 (48.1)	16 (45.7)
Grade III	2 (7.5)	3 (9.6)
Metastasis at diagnosis, n (%)			>0.05
Yes	2 (7.5)	4 (11.4)
No	25 (92.5)	31 (88.6)
Tumor subsite, n (%)			>0.05
Oral tongue	10 (37.0)	13 (37.1)
Buccal mucosa	2 (7.4)	2 (5.7)
Larynx	7 (26.0)	9 (25.7)
Oral cavity	4 (14.8)	5 (14.2)
Floor of mouth	4 (14.8)	6 (17.1)

[i] HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma.

Data processing

For constructing the WGCNA, the gene expression profiles of 54 HPV-negative HNSCC tissues (ISH and P16 negative stain), along with their corresponding clinical characteristics (Table SI), were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). After normalization, outliers were assessed using the hierarchical cluster algorithm with a cut-off of 140. A fragments per kilobase of transcript per million mapped reads (FPKM) value of <0.5 was used as the threshold for removing the genes with low expression. This was necessary as a gene must be expressed at a minimal level before it is likely to be translated into a protein or be considered biologically important (13). Consequently, 15,703 genes expressed in the 54 HPV-negative HNSCC tissues, along with the corresponding clinical characteristics data, were used in the WGCNA. For DEGs analysis, the gene expression prolife in TCGA and Gene Expression Omnibus datasets [GSE117973 (14), GSE85446 (15) and GSE112026 (16); https://www.ncbi.nlm.nih.gov/gds] containing 103 HPV-positive HNSCC and 165 HPV-negative HNSCC was downloaded. After merge and batch normalization using the sva package (version: 3.11; Bioconductor; https://www.bioconductor.org/packages/release/bioc/html/sva.htl), the merge gene expression profile were used for DEGs analysis. |Log fold-change (FC)| >0.5 and adjusted P<0.05 were set as cut-off values to consider genes differentially expressed between HPV-positive and HPV-negative HNSCC.

WGCNA

The WGCNA R package (version 1.69; https://cran.r-project.org/web/packages/WGCNA/index.html) was used to conduct the WGCNA. All gene pairs were analyzed using Pearson's correlation analysis, the results of which were used to construct a matrix of similarity. Subsequently, to produce a scale-free co-expression network, the matrix of similarity was constructed using a soft power of β=7. The adjacency matrix was then translated into a topological overlap matrix (TOM). Furthermore, median linkage hierarchical clustering was analyzed using the TOM-based dissimilarity measure with a minimum size of 50.

Identification of clinically significant modules

Following WGCNA, the Pearson's correlation analyses between various module eigengenes and the clinical characteristics of patients were assessed. These traits included node (N) stage, tumor (T) stage, tumor grade, recurrence, perineural invasion, tobacco use history and alcohol consumption history. Modules were considered clinically significant if they were correlated with two clinical characteristics (r>0.3 and P<0.05). In the clinically significant modules, gene significance (GS) was quantified using associations between the individual genes and the clinical characteristics of interest, along with the module membership (MM), which was itself determined using the correlation between the module eigengenes and the gene expression profiles. If P<0.05 for the correlation between GS and MM in the clinically significant modules, the modules were subjected to further analysis.

Construction of a protein-protein interaction (PPI) network

The online Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://string-db.org) was used to construct the PPI network. Nodes without connections were removed. The minimum required interaction score was set at 0.7 (high confidence). Node and edge information was exported into a text file, which was subsequently imported into Cytoscape software (version 6.1; http://cytoscape.org/). In the PPI network, nodes represented the protein information encoded by the gene of interest. Lines indicated the interactions between two proteins. The degree score was calculated using the Cytohub plug-in (version 1.0; http://github.com/cytoscape/appstore) and was set as the criterion used to calculate the interaction between a single protein and other proteins. Finally, proteins encoded by genes with a degree score in the top 10% were visualized using Cytoscape, and the genes coding for proteins with a degree score in the top 1% were considered candidate hub genes.

Pearson's correlation analysis between the candidate hub genes and tumor grade

Pearson's correlation analysis between the expression of candidate hub genes and tumor grade was performed using SPSS version 20.0 (IBM Corp.). P<0.05 was set as the cut-off indicating a significant correlation. The genes significantly correlated with tumor grade were considered hub genes.

Immunohistochemical (IHC) staining

All 27 HPV-negative and 35 HPV-positive HNSCC tissues were fixed in 4% paraformaldehyde (Boster Biological Technology) at room temperature. Subsequently, they were embedded into paraffin and cut into 4-µm-thick sections. After heating at 65°C, the sections were deparaffinized using xylene and rehydrated in a graded ethanol series (100, 80, 60 and 40%). After conducting antigen retrieval using sodium citrate (0.01 mol/l; Boster Biological Technology), endogenous peroxidase activity was blocked using 3% H2O2. Subsequently, 5% bovine serum albumin (Wuhan Servicebio Technology Co., Ltd.) was added at room temperature for 1 h. The sections were then incubated for 16 h at 4°C with the following primary antibodies: p16INK4a (1:200; cat. no. A11651; ABclonal Biotech Co., Ltd.), Serrate RNA effector molecule (SRRT; 1:200; cat. no. A8219; ABclonal Biotech Co., Ltd.), checkpoint kinase 2 (CHEK2; 1:200; cat. no. A19543; ABclonal Biotech Co., Ltd.), small nuclear ribonucleoprotein polypeptide E (SNRPE; 1:200; cat. no. A5488; ABclonal Biotech Co., Ltd.), proteasome 26S subunit ATPase 2 (PSMC2; 1:400; cat. no. 14905-1-AP; ProteinTech Group, Inc.), origin recognition complex subunit 5 (ORC5; 1:400; cat no. 11542-1-AP, ProteinTech Group, Inc.), S100 calcium binding protein A7 (S100A7; 1:200; cat no. 13061-1-AP, ProteinTech Group Inc.) and keratinocyte differentiation associated protein (KRTDAP; 1:1,000; cat. no. ab204583; Abcam). Subsequently, the sections were incubated for 2 h at room temperature with anti-mouse (cat. no. BM3895) and anti-rabbit (cat. no. BM3894) horseradish peroxidase-conjugated goat secondary antibodies (both 1:500; Boster Biological Technology). Samples were then stained with 3,3′-diaminobenzidine and hematoxylin at room temperature for 1 min, after which a light orthophoto microscope (magnification, ×200 and ×400) was used to obtain images. Finally, the protein expression levels of the target genes were evaluated based on the sum of the intensity score (0, no staining; 1, weakly positive; 2, moderately positive; and 3, strongly positive) and the score for the proportion of positive cells (0, <1; 1, 1–33; 2, 34–66; and 3, 67–100%) using Image-Pro Plus software (version 6.0; Media Cybernetics, Inc.). Total scores of 0–2, 3–4 and 5–6 indicated low, moderate and high expression, respectively. The difference in the protein expression of tissues was determined using Fisher's exact test combined with Bonferroni's correction test based on their expression levels. P<0.05 was considered to indicate a statistically significant difference.

Verification of the diagnostic value of hub genes

The diagnostic value of hub genes for HPV-negative and HPV-positive HNSCC was analyzed using a receiver operating characteristic (ROC) curve analysis using the IHC-based protein levels. The protein expression scores of the HNSCC and normal tissues were imported into SPSS version 20.0 (IBM Corp.), after which ROC curve analysis was performed. An area under the curve (AUC) value of >0.7 was considered to indicate a high diagnostic value.

Results

WGCNA

The gene expression profiles of 54 HPV-negative HNSCC tissues (with negative ISH and P16 staining) were downloaded from the TCGA database, as well as corresponding clinical characteristic data. The sample dendrogram revealed that there were no outliers; similarly, the trait heatmap showed that clinical data of most patients in TCGA database were completely documented (Fig. 1). Therefore, after removing the genes with low expression, the gene expression profiles of all 54 HPV-negative HNSCC tissues were used in the WGCNA. When the soft power β was set as 7, the scale independence of the topology network reached >0.85 (Fig. 2A) with the mean connectivity close to 0 (Fig. 2B). Therefore, a soft power of β=7 was selected as the soft threshold for performing subsequent analyses. The results revealed that when β=7, the topological overlap matrix was able to meet the scale-free topology criterion with R2=0.87 (Fig. 2C). The results identified 27 gene co-expression modules (paleturquoise, cyan, salmon, sienna3, violet, royalblue, saddlebrown, darkmagenta, brown, black, purple, lightyellow, darkolivegreen, grey60, yellow, darkgreen, skyblue, lightcyan, lightgreen, tan, darkturquoise, darkorange, blue, green, steelblue, darkgrey and midnightblue), while the genes that were not co-expressed were clustered in the grey module (Fig. 2D). The 27 gene co-expression modules were then used for further analysis.

Figure 1. Sample dendrogram of 54 human papillomavirus-negative head and neck squamous cell carcinoma tissues and the corresponding clinical characteristics of patients, including N stage, T stage, tumor grade, recurrence, perineural invasion, tobacco use history and alcohol consumption history. Red, white and grey indicate high, low and missing values, respectively. N, node; T, tumor.

Figure 2. Weighted gene co-expression network analysis. (A) Analysis of the scale-free fit index for various β soft-thresholding powers. (B) Analysis of the mean connectivity for various soft-thresholding powers. (C) Scale-free topology when β=7. (D) Dendrogram of all genes clustered based on a dissimilarity measure (1-topological overlap matrix). There were 27 gene co-expression modules in total (paleturquoise, cyan, salmon, sienna3, violet, royalblue, saddlebrown, darkmagenta, brown, black, purple, lightyellow, darkolivegreen, grey60, yellow, darkgreen, skyblue, lightcyan, lightgreen, tan, darkturquoise, darkorange, blue, green, steelblue, darkgrey and midnightblue). Genes without co-expression were clustered into the grey module.

Identification of clinically significant modules

By analyzing the correlation between the gene co-expression modules and clinical characteristics, it was determined that four gene co-expression modules were simultaneously and significantly correlated with two clinical characteristics in patients with HPV-negative HNSCC. The gene of the grey60 module was positively correlated with T stage (R=0.38; P=0.004) and tumor grade (R=0.32; P=0.03). Genes of the lightyellow module were positively correlated with tumor grade (R=0.31; P=0.04) and negatively correlated with perineural invasion (R=−0.4; P=0.002). Additionally, genes in the purple module were positively correlated with T stage (R=0.32; P=0.02) and tumor grade (R=0.31; P=0.03), while genes of the salmon module were negatively correlated with T stage (R=−0.31; P=0.02) and tumor stage (R=−0.39; P=0.003) (Fig. 3). Correlations between GS and MM were subsequently calculated in the four aforementioned modules. The MM of the grey60 module was correlated with the GS for tumor grade (correlation=0.29; P<0.05) and T stage (correlation=0.31; P<0.05) (Fig. 4A). The MM of the lightyellow module was correlated with the GS for tumor grade (correlation=0.35; P<0.05) and perineural invasion (correlation=0.7; P<0.05) (Fig. 4B). The MM of the purple module was correlated with the GS for tumor grade (correlation=0.3; P<0.05) and T stage (correlation=0.45; P<0.05) (Fig. 4C), while the MM of the salmon module was correlated with the GS for tumor grade (correlation=0.35; P<0.05) and T stage (correlation=0.32; P<0.05) (Fig. 4D). Therefore, these four modules (grey60, lightyellow, purple and salmon modules) were set as key gene modules and were analyzed further.

Figure 3. Identification of clinically significant modules based on the correlation of module eigengenes with clinical characteristics. Each cell contains the corresponding correlation coefficient and P-value (in brackets). The colors indicate the strength of correlation according to the color legend. ME, module; T, tumor; N, node.

Figure 4. Analysis of GS and MM in the clinically significant modules. (A) MM of the grey60 module was correlated with GS for tumor grade and T stage. (B) MM in the lightyellow module was correlated with GS for tumor grade and perineural invasion. (C) MM of the purple module was correlated with GS for tumor grade and T stage. (D) MM of the salmon module was correlated with GS for tumor grade and T stage. GS, gene significance; MM, module membership; cor, correlation; T, tumor; N, node.

Construction of the PPI network and selection of candidate hub genes

The lightyellow, grey60, purple and salmon modules included 201, 997, 2,538 and 389 genes, respectively. After importing these genes into the STRING database and removing the nodes without connections, genes with a degree score in the top 10% in the grey60 (Fig. 5A), lightyellow (Fig. 5B), purple (Fig. 5C) and salmon (Fig. 5D) modules were visualized using Cytoscape. The genes with a degree score in the top 1% in each module are presented in Table II. The resultant 17 genes were considered candidate hub genes and may serve key roles in the development of HPV-negative HNSCC.

Figure 5. PPI network. Genes with a degree score in the top 10% in the (A) grey60, (B) light yellow, (C) purple and (D) salmon modules were used to construct the PPI network. Nodes and edges indicate the genes and their interactions, respectively. PPI, protein-protein interaction.

Table II. Genes with a degree score in the top 1% in each module.

Table II.

Genes with a degree score in the top 1% in each module.

Module	Gene	Degree score
Grey60 module	SNRPG	20
	SNRPF	18
	SNRPE	16
	SNRPD1	16
	NOP58	16
Salmon module	S100A7	8
	KRTDAP	7
Purple module	DDX39B	5
	CHEK2	4
	KAT2A	4
Lightyellow module	ORC5	4
	PSMC2	4
	SRRT	4
	BUD31	4
	POLR2J	4
	ARMC10	4
	RBM48	4

[i] SRRT, serrate RNA effector molecule; CHEK2, checkpoint kinase 2; SNRPE, small nuclear ribonucleoprotein polypeptide E; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; S100A7, S100 calcium binding protein A7; KRTDAP, keratinocyte differentiation associated protein.

Pearson's correlation analysis between the candidate hub genes and tumor grade

As the four gene modules were all correlated with tumor grade, the correlation between candidate hub gene expression and the tumor grade of HPV-negative HNSCC was further analyzed using Pearson's correlation analysis. The results revealed that SRRT (R=0.30; P<0.05), CHEK2 (R=0.29; P<0.05), SNRPE (R=0.31; P<0.05), PSMC2 (R=0.31; P<0.05) and ORC5 (R=0.32; P<0.05) were positively correlated with the tumor grade of HPV-negative HNSCC tissue, while S100A7 (R=−0.38; P<0.05) and KRTDAP (R=−0.29; P<0.05) were negatively correlated with tumor grade (Fig. 6). Therefore, SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP were identified as hub genes in HPV-negative HNSCC.

Figure 6. Pearson's correlation analysis between hub gene expression and tumor grade in human papillomavirus-negative head and neck squamous cell carcinoma. SRRT, serrate RNA effector molecule; CHEK2, checkpoint kinase 2; SNRPE, small nuclear ribonucleoprotein polypeptide E; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; S100A7, S100 calcium binding protein A7; KRTDAP, keratinocyte differentiation associated protein.

Identification of the expression levels of hub genes in HPV-negative and HPV-positive HNSCC

To analyze the expression levels of hub genes in HPV-negative and HPV-positive HNSCC, the merged gene data expression profile of 103 HPV-positive HNSCC and 165 HPV-negative HNSCC samples in TCGA and GEO databases (GSE117973, GSE85446 and GSE112026) was used to perform DEGs analysis. The results revealed that there were 27 upregulated genes and 25 downregulated genes in HPV-negative HNSCC compared with in HPV-positive HNSCC (Fig. 7A; Table III). Through an intersection analysis, it was found that the hub genes ORC5 and PSMC2 were upregulated in HPV-negative HNSCC tissues compared with in HPV-positive HNSCC tissues, while KRTDAP was downregulated (Fig. 7B).

Figure 7. Differential expression gene analysis for HPV-negative and HPV-positive HNSCC. (A) Volcano plot of the differentially expressed genes between HPV-negative and HPV-positive HNSCC. (B) Intersection analysis for the differentially expressed genes (HPV-negative HNSCC versus HPV-positive HNSCC) and hub genes identified from WGCNA. HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma; WGCNA, weighted gene co-expression network analysis; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; KRTDAP, keratinocyte differentiation associated protein.

Table III. Detailed information of genes differently expressed in HPV-negative and HPV-positive head and neck squamous cell carcinoma according to the different expression analysis.

Table III.

Detailed information of genes differently expressed in HPV-negative and HPV-positive head and neck squamous cell carcinoma according to the different expression analysis.

Gene	Log (fold-change)
DSG1	0.814689016
SPRR2G	0.807843405
KLK5	0.78303004
KRT1	0.735022158
PSMC2	0.729090901
MMP10	0.701782465
ORC5	0.655262803
SPRR2B	0.650839232
LCE3E	0.646742549
DSC1	0.622678824
KLK7	0.622523674
DEFB103B	0.602302594
PTHLH	0.598655999
MMP13	0.588604647
MMP3	0.584023444
FAM25A	0.576520859
ASPRV1	0.566257594
LCE3D	0.564179626
CCNA1	0.561757208
KLK8	0.558113758
CPA4	0.555893249
IGFL1	0.551152103
CDSN	0.532585037
KLK6	0.52189567
CDA	0.512697702
KRT75	0.507135684
WFDC12	0.506385839
CYP4X1	−0.50678768
SYNGR3	−0.514548202
RIBC2	−0.514701295
RNF212	−0.517145557
ABCA3	−0.528349222
TMPRSS2	−0.530093548
TMSB15A	−0.53215901
KRTDAP	−0.54853466
MEI1	−0.551064139
GABRP	−0.554796624
CDKN2C	−0.565098291
YBX2	−0.569455778
TCP11	−0.570473276
KRT19	−0.578420257
SMC1B	−0.580791027
PODXL2	−0.59651466
PLAC8	−0.620386125
FAM3B	−0.630548652
STAG3	−0.665598925
TAF7L	−0.672735588
ZNF541	−0.677964463
KCNS1	−0.683860665
SYCP2	−0.735795983
NEFH	−0.757761211
CDKN2A	−0.768176467

[i] HPV, human papillomavirus; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; KRTDAP, keratinocyte differentiation associated protein.

Verification of the protein expression levels of hub genes in HPV-negative and HPV-positive HNSCC tissues

To verify the bioinformatics results, IHC staining of 27 HPV-negative and 35 HPV-positive HNSCC tissues was performed, alongside their corresponding adjacent tissues. The results revealed that the protein expression levels of SRRT, CHEK2 and SNRPE were increased in HPV-negative HNSCC and HPV-positive HNSCC tissues compared with in their adjacent counterparts; however, there was no significant difference between HPV-negative and HPV-positive HNSCC tissues. S100A7 expression was not significantly different among HPV-negative HNSCC, HPV-positive HNSCC and adjacent tissues. Furthermore, it was demonstrated that PSMC2 and ORC5 expression was significantly increased, while KRTDAP expression was significantly decreased in HPV-negative HNSCC tissues compared with in adjacent and HPV-positive HNSCC tissues (Table IV; Figs. 8 and S2).

Figure 8. Immunohistochemistry was used to detect the protein expression levels of SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP in HPV-negative HNSCC, HPV-positive HNSCC and adjacent tissues (magnification, ×200 and ×400). Black scale bar, 100 µm; white scale bar, 50 µm. SRRT, serrate RNA effector molecule; CHEK2, checkpoint kinase 2; SNRPE, small nuclear ribonucleoprotein polypeptide E; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; S100A7, S100 calcium binding protein A7; KRTDAP, keratinocyte differentiation associated protein; HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma; N, normal tissues; T, tumor tissues.

Table IV. Detailed immunohistochemical scoring of SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP expression in HPV− HNSCC (n=27), HPV+ HNSCC (n=35) and corresponding adjacent tissues.

Table IV.

Detailed immunohistochemical scoring of SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP expression in HPV− HNSCC (n=27), HPV+ HNSCC (n=35) and corresponding adjacent tissues.

Gene	Tissues	Low expression (0–2), n	Medium expression (3–4), n	High expression (5–6), n
SRRT	Tumor (HPV+)a	5	4	26
	Adjacent (HPV+)	27	6	2
	Tumor (HPV−)b	2	3	22
	Adjacent (HPV−)	18	6	3
CHEK2	Tumor (HPV+)a	3	3	29
	Adjacent (HPV+)	28	4	3
	Tumor (HPV−)b	5	4	18
	Adjacent (HPV−)	18	7	2
PSMC2	Tumor (HPV+)	15	13	7
	Adjacent (HPV+)	22	7	6
	Tumor (HPV−)b,c	2	3	22
	Adjacent (HPV−)	20	5	2
OCR5	Tumor (HPV+)	19	10	6
	Adjacent (HPV+)	24	9	2
	Tumor (HPV−)b,c	2	7	18
	Adjacent (HPV−)	19	5	3
S100A7	Tumor (HPV+)	6	5	24
	Adjacent (HPV+)	4	4	27
	Tumor (HPV−)	2	4	21
	Adjacent (HPV−)	4	5	18
SNRPE	Tumor (HPV+)a	3	4	28
	Adjacent (HPV+)	26	3	6
	Tumor (HPV−)b	4	2	21
	Adjacent (HPV−)	17	6	4
KRTDAP	Tumor (HPV+)	9	6	20
	Adjacent (HPV+)	8	8	19
	Tumor (HPV−)b,c	20	4	3
	Adjacent (HPV−)	4	7	16

a P<0.05 vs. HPV⁺ adjacent tissues

b P<0.05 vs. HPV⁻ adjacent tissues

c P<0.05 vs. HPV+ HNSCC. SRRT, serrate RNA effector molecule; CHEK2, checkpoint kinase 2; SNRPE, small nuclear ribonucleoprotein polypeptide E; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; S100A7, S100 calcium binding protein A7; KRTDAP, keratinocyte differentiation associated protein; HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma.

Verification of the diagnostic value of hub genes for HPV-negative and HPV-positive HNSCC

Using protein level data obtained by IHC staining, ROC curve analysis was performed. The results revealed that SRRT, CHEK2, SNRPE, PSMC2, ORC5 and KRTDAP exhibited high diagnostic value for distinguishing between HPV-negative HNSCC and adjacent tissues, while S100A7 did not exert this effect (Fig. 9A). However, only SRRT, CHEK2 and SNRPE had high diagnostic value for distinguishing between HPV-positive HNSCC and adjacent tissues (Fig. 9B). Thus, SRRT, CHEK2 and SNRPE may be used as common biomarkers for both HPV-negative and HPV-positive HNSCC, while PSMC2, ORC5 and KRTDAP may be specific biomarkers for HPV-negative HNSCC only.

Figure 9. Receiver operating characteristic curve analysis was used to determine the diagnostic value of SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP for (A) HPV-negative and (B) HPV-positive HNSCC based on protein level data obtained by immunohistochemical staining. AUC >0.7 indicated genes with high diagnostic value for distinguishing between HNSCC and adjacent tissues. SRRT, serrate RNA effector molecule; CHEK2, checkpoint kinase 2; SNRPE, small nuclear ribonucleoprotein polypeptide E; PSMC2, proteasome 26S subunit ATPase 2; ORC5, origin recognition complex subunit 5; S100A7, S100 calcium binding protein A7; KRTDAP, keratinocyte differentiation associated protein; HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma; AUC, area under the curve.

Discussion

There are numerous risk factors involved in the pathogenesis of HNSCC, including HPV infection. Although the majority of HNSCC cases involve HPV infection, certain patients develop HNSCC without it (17,18). The molecular mechanisms underlying these two subtypes are different (19). Furthermore, HPV-negative HNSCC is more aggressive than HPV-positive HNSCC (20,21). Therefore, the identification of novel and specific biomarkers for HPV-negative HNSCC may improve the understanding of the specific molecular mechanism associated with HPV-negative HNSCC, which may improve the diagnosis and treatment of patients.

The current study identified four gene co-expression modules associated with clinical characteristics of patients via WGCNA. The results revealed 17 genes in these gene co-expression modules that had a high degree score (top 1% of genes in each module) in the PPI network. Additionally, seven of these genes (SRRT, CHEK2, SNRPE, PSMC2, ORC5, S100A7 and KRTDAP) were correlated with tumor grade in HPV-negative HNSCC. DEG analysis demonstrated that PSMC2 and ORC5 expression was higher in HPV-negative HNSCC compared with in HPV-positive HNSCC, while KRTDAP expression was lower. Furthermore, using IHC staining and ROC curve analysis, it was revealed that SRRT, CHEK2, SNRPE, PSMC2, ORC5 and KRTDAP were differentially expressed in HPV-negative HNSCC tissues compared with in adjacent tissues and had high diagnostic values for distinguishing between HPV-negative HNSCC tissues and adjacent tissues. However, after performing additional IHC staining and ROC curve analysis, it was demonstrated that SRRT, CHEK2 and SNRPE were differentially expressed in HPV-positive HNSCC tissues compared with in adjacent tissues, exhibiting a high diagnostic value for distinguishing between HPV-positive HNSCC and adjacent tissues. Thus, while SRRT, CHEK2 and SNRPE may be used as common biomarkers for both HPV-negative and HPV-positive HNSCC, PSMC2, ORC5 and KRTDAP may be specific to HPV-negative HNSCC only.

SRRT (also named Ars2) was first isolated from a hamster cell line and serves a key role in sodium arsenite resistance (22). Furthermore, it has been demonstrated that SRRT is involved in the biosynthesis of certain microRNAs (23). However, research on the role of SRRT in the development of cancer is conflicting. It has been demonstrated that SRRT serves as an oncogene in glioblastoma and promotes the proliferation and migration of LN-229 and U87 glioblastoma cells via the MAPK/ERK signaling pathway (24). Additionally, SRRT expression is upregulated in cholangiocarcinoma and associated with poor outcomes (25). However, pediatric acute lymphoblastic leukemia involving low SRRT expression has poor chemotherapy outcomes (26). CHEK2 was first reported as an important breast cancer susceptibility gene; it serves a key role in regulating cell apoptosis, cell cycle and DNA repair (27). It increases the stability of P53 and induces cell cycle arrest in the G1 phase (28). However, CHEK2 mutations are very common in several types of cancer, including breast cancer (29), colorectal cancer (30) and oral squamous cell carcinoma (31). Additionally, CHEK2 promotes the progression of various types of cancer, including hepatocellular carcinoma (32) and colorectal cancer (33). SNRPE is a member of a large family of polypeptides that are conserved in eukaryotes and archaebacteria (34). Previous studies have revealed that it is involved in RNA processing and mRNA degradation (35,36). Furthermore, a previous study has revealed that SNRPE is highly expressed in high-grade prostate cancer, promoting cell proliferation (37). In the present study, it was demonstrated that SRRT, CHEK2 and SNRPE were in clinically significant gene modules, had high degree scores in the PPI network, were associated with tumor grade and were highly expressed in HPV-negative and HPV-positive HNSCC tissues. Furthermore, it was revealed that SRRT, CHEK2 and SNRPE had a high diagnostic value in distinguishing HNSCC and adjacent tissues. Therefore, to the best of our knowledge, the present study was the first to demonstrate that SRRT, CHEK2 and SNRPE may be used as common biomarkers for both HPV-negative and HPV-positive HNSCC. Targeting SRRT, CHEK2 and SNRPE may aid the therapy of patients with HPV-negative and HPV-positive HNSCC. However, the mechanism of SRRT, CHEK2 and SNRPE in the development of HNSCC requires further study.

PSMC2 is a novel gene on chromosome 7q22.1-q22.3 that encodes a member of the 19S proteasome (38). Recently, PSMC2 was revealed as an oncogene in several types of cancer, including pancreatic cancer (39) and osteosarcoma (40). Similarly, high PSMC2 expression predicts a poor prognosis in patients with osteosarcoma (41). ORC5 is a key member of the origin recognition complex that binds to replication origins, initiating replication itself (42). Additionally, it induces large-scale chromatin decondensation and regulates the cell cycle and proliferation (43). KRTDAP acts as a soluble regulator of keratinocyte differentiation and serves a key role in embryonic skin morphogenesis (44). Research has revealed that KRTDAP mRNA is upregulated in HPV-positive HNSCC tissues compared with in adjacent tissues (45). In the present study, it was demonstrated that PSMC2 and ORC5 expression was significantly increased in HPV-negative HNSCC tissues compared with in corresponding adjacent and HPV-positive HNSCC tissues. Additionally, KRTDAP was significantly decreased in HPV-negative HNSCC tissues compared with in corresponding adjacent and HPV-positive HNSCC tissues. Similarly, it was revealed that PSMC2, ORC5 and KRTDAP had a high diagnostic value in distinguishing HPV-negative HNSCC and adjacent tissues. However, it could not distinguish between HPV-positive HNSCC and adjacent tissues. Therefore, it was hypothesized that PSMC2, ORC5 and KRTDAP may be used as specific biomarkers for patients with HPV-negative HNSCC only. It is well known that HPV-negative HNSCC is more aggressive than HPV-positive HNSCC, so the downregulation of KRTDAP and upregulation of PSMC2 and ORC5 may be key for HPV-negative HNSCC cells to become aggressive. Targeting ORC5, PSMC and KRTDAP may therefore help the individualized treatment of patients with HPV-negative HNSCC. However, additional experiments should be performed to verify this hypothesis.

In conclusion, based on WGCNA, PPI network, Pearson's correlation, DEG, IHC and ROC curve analyses, SRRT, CHEK2 and SNRPE were revealed to be common biomarkers for both HPV-negative and HPV-positive HNSCC, while PSMC2, ORC5 and KRTDAP may be potential specific biomarkers for HPV-negative HNSCC. The results may improve the understanding of HNSCC, as well as providing potential biomarkers and targets for the diagnosis and treatment of patients with this disease.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The current study was supported by Guizhou Big Data Health Management Industry Technology Innovation Strategic Alliance [grant no. Qianketong (2016)193].

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request, and the GSE117973, GSE85446 and GSE112026 datasets analyzed during the current study are available in the Gene Expression Omnibus repository (https://www.ncbi.nlm.nih.gov/gds).

Authors' contributions

YS, YC and ZZ were responsible for the experimental design, analysis and interpretation of data. DR, YY and BW performed the experiments. YS and YC were responsible for confirming the authenticity of the data. All authors took part in drafting the article and revising it critically for important intellectual content, agreed to be accountable for all aspects of the work, and read and approved the final manuscript.

Ethics approval and consent to participate

The current study was approved by the Human Trait Ethics Committee of Guizhou Medical University (Guiyang, China) and performed in accordance with the Declaration of Helsinki. All patients from which samples were obtained provided written informed consent.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References