A HPV16-related prognostic indicator for head and neck squamous cell carcinoma.

BACKGROUND: The human papillomavirus (HPV) is emerging as an important risk factor in head and neck squamous cell carcinoma (HNSCC) patients. This has been observed particularly in the case of HPV16. The HPV16+ HNSCC subtype has distinct pathological, clinical, molecular, and prognostic characteristics. This study aimed to identify potential microRNAs (miRNAs) and their roles in HPV16+ HNSCC progression.
METHOD: miRNA, mRNA and the clinical data of 519 HNSCC and 44 HNSCC-negative samples were obtained from The Cancer Genome Atlas (TCGA) database. Differentially expressed miRNAs (DEMs) in HPV16-related HNSCC tissues with prognostic value were selected. DEM levels were assessed based on clinicopathological parameters and overall survival (OS). Target genes were also predicted and functional analysis based on Gene Set Enrichment Analysis (GSEA) were then performed.
RESULTS: In HPV16+ HNSCC tissues, miR-99a-3p and miR-4746-5p were significantly upregulated. In contrast, miR-411-5p was shown to be downregulated. miR-99a-3phighmiR-411-5plowmiR-4746-5phigh expression could estimate improved OS and low frequent perineural invasion (PNI). Predicted target genes were enriched in cell growth, neuroepithelial cell differentiation, MAPK and FoxO signaling pathways. Epithelial mesenchymal transition (EMT) gene set and invasion related genes were downregulated in miR-99a-3phighmiR-411-5plowmiR-4746-5phigh HNSCC patients.
CONCLUSION: miR-99a-3p, miR-411-5p and miR-4746-5p might participate in HPV16+ HNSCC progression through EMT related pathways and affect prognosis. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

The human papillomavirus (HPV) is a major cause for infectious agent-related cancers, including cervical cancer and head and head and neck squamous cell carcinoma (HNSCC) (1). It contributes to an estimated 600,000 or more incident cancer cases and 250,000 deaths each year (2,3). In addition, an increasing number of studies have demonstrated superior prognosis and different genomic and phenotypic characteristics of patients with HPV-associated HNSCC (4-7).

MicroRNAs (miRNAs) play pivotal roles in post-transcriptional regulation of gene expression in HPV-related cancers (8,9). High-risk HPV16 type genes (E5, E6, E7) have been reported to interact with miRNAs in tumor progression (10-13). HPV+ HNSCC showed distinct miRNA profile. These miRNAs may potentially target common HPV-mediated pathways such as E6/p53 and E7/pRb (14,15). However, miRNA function on the prognosis of HPV16+ HNSCC was unclear.

This study aimed to identify potential miRNAs and their roles in HPV16+ HNSCC progression, based on the The Cancer Genome Atlas (TCGA) database. Potential target genes were predicted, and possible functions were also analyzed. We present the following article in accordance with the REMARK reporting checklist (available at http://dx.doi.org/10.21037/atm-20-6338).

Methods

Data processing

All clinical and sequencing data was obtained from The TCGA database (https://cancergenome.nih.gov/). Based on the inclusion criteria of cases with complete miRNA expression, mRNA expression and corresponding clinical data, 519 HNSCC patients and 44 corresponding HNSCC-negative control samples were included in this study. The raw data was converted into a data matrix in R, and the identification of differentially expressed miRNAs (DEMs) and mRNAs (DEGs) was performed in HNSCC and normal samples through the R package, Limma. In addition, DEMs between HNSCC patients with different HPV 16 statuses were also selected. All DEMs or DEGs meeting the cutoff criteria |log2(Fold Change)|>1 and P<0.05 were considered significant. Univariate Cox regression analysis was performed to select significant prognostic miRNAs. Common DEMs were identified using the plug-in jvenn (http://jvenn.toulouse.inra.fr/app/index.html) (16), and the final prognostic miRNAs were determined by multivariate Cox regression analysis. In addition, a risk model based on 3 miRNAs was built with a risk . All patients were divided into high- or low-risk groups by the rank of risk value.

Correlation of DEMs with clinical and pathological characteristics of HNSCC patients

The expression levels of selected DEMs with clinical and pathological parameters (gender, age, alcohol and smoking history, cancer grade, cancer stage, T-stage, N stage, HPV16 status, perineurial invasion (PNI), angiolymphatic invasion, disease surgical margin status) were compared using the unpaired t test or one-way ANOVA.

Association of DEMs and patient prognosis

The expression levels of selected DEMs underwent log2 transformation before analysis. Each patient was classified as having either positive or negative expression according to the corresponding median miRNA expression. The Kaplan-Meier approach with log-rank analysis using the R package (survival) was adopted to determine the role of candidate prognostic miRNAs in HNSCC. In addition, overall survival (OS) analysis for HNSCC patients with high or low risk, based on the calculated risk score was also performed to estimate the prognostic value of all selected DEMs.

Target gene prediction

Predicted target genes of the selected DEMs were verified using the miRDB database (17).

Gene Ontology (GO) functional and pathway enrichment analyses and Gene Set Enrichment Analysis (GSEA)

The GO functions with respect to the biological process (BP), molecular function (MF), and cellular component (CC) categories and KEGG pathway enrichment for the predicted target genes were identified through the Metascape database with a P<0.05 cutoff.

GSEA for DEGs in HPV16+ HNSCC patients [|log2(Fold Change)|>1 and P<0.05] was performed using R package clusterProfiler. Hallmark gene sets were selected for functional analysis.

Statistics analysis

Statistical analyses were performed using the R program. Analytical data is expressed as mean ± standard deviation (SD). Unpaired t-tests or one-way ANOVA were performed to assess the relationships between clinical and pathological characteristics, along with miRNA expression levels. Analysis of parameters with 2 groups was performed using an unpaired t-test, while analysis of parameters with 3 or more groups was performed using one-way ANOVA. Pearson correlation analysis was performed for each pair of miRNAs expressed in HNSCC samples. Univariate and multivariate Cox proportional hazard regression analyses and Kaplan-Meier survival analysis were conducted to compare the prognostic values of high- and low-level miRNAs. Differences with a P value <0.05 were considered significant.

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Determination of prognostic DEMs in HPV+ HNSCC

lists the expression of miR-99a-3p, miR-411-5p and miR-4746-5p in HNSCC samples based on clinicopathological parameters. Based on the cutoff criteria (P<0.05 and |log2FC|>1.0), 183 DEMs were identified between cancer and normal samples. An additional 59 DEMs were found between HNSCC patients with varied HPV16 statuses (). Furthermore, 195 potential prognostic miRNAs were selected by univariate analysis with the cutoff criterion of P<0.05. As shown in , 6 common miRNAs were selected (miR-99a-3p, miR-411-5p, miR-432-5p, miR-495-3p, miR-758-5p, miR-4746-5p). Subsequently, final prognostic miRNAs (miR-99a-3p, miR-411-5p, miR-4746-5p) were confirmed by multivariate Cox regression analysis with the cutoff P<0.05 ().

Table 1

Expression of miR-99a-3p, miR-411-5p and miR-4746-5p in HNSCC samples, based on clinicopathological parameters

Parameter	Group	Cases	miR-99a-3p				miR-411-5p				miR-4746-5p
			T/F	Expression (log2 RPM +0.01, ± SD)	P value		T/F	Expression (log2 RPM +0.01, ± SD)	P value		T/F	Expression (log2 RPM +0.01, ± SD)	P value
Age	≥60	287	3.205	1.154, 2.307	0.001		1.448	2.064, 1.523	0.148		−0.239	2.160, 1.196	0.811
	<60	231		0.476, 2.511				1.896, 1.052				2.186, 1.333
Gender	Male	379	0.196	0.472, 2.347	0.845		−0.344	2.588, 1.394	0.731		3.328	2.273, 1.171	0.001
	Female	140		0.431, 4.240				2.634, 1.084				1.867, 1.508
Alcohol	Yes	349	−0.506	0.411	0.613		1.250	2.652, 1.338	0.212		1.250	2.652, 1.338	0.212
	No	161		0.516				2.485, 1.423				2.485, 1.423
Smoke	Non-smoker	121	−0.848	0.293, 2.502	0.398		-2.137	2.376, 1.346	0.034		-0.354	2.136, 1.094	0.724
	Smoker	386		0.508, 2.210				2.676, 1.368				2.178, 1.240
HPV p16	Positive	40	2.130	1.263, 2.435	0.037		-4.110	1.614, 1.429	<0.001		4.155	3.035, 1.228	<0.001
	Negative	73		0.306, 1.977				2.739, 1.320				2.011, 1.297
T stage	T1	38	4.262	0.865, 1.964	0.005		0.104	2.709, 1.588	0.958		0.202	2.268, 1.155	0.895
	T2	147		0.878, 2.040				2.584, 1.415				2.193, 1.419
	T3	139		−0.008, 2.477				2.640, 1.348				2.124, 1.004
	T4	183		0.430, 2.170				2.638, 1.252				2.136, 1.178
N stage	N0	243	1.204	0.591, 2.064	0.308		1.354	2.635, 1.260	0.256		1.479	2.098, 1.248	0.220
	N1	83		0.116, 2.304				2.831, 1.457				2.035, 1.143
	N2	165		0.438, 2.454				2.516, 1.427				2.315, 1.161
	N3	9		1.136, 2.001				3.090, 1.303				2.002, 1.253
M status	M0	493	−0.056	0.479, 2.225	0.957		0.500	2.629, 1.365	0.638		0.388	2.153, 1.211	0.714
	M1	6		0.513, 1.447				2.345, 1.384				1.960, 1.207
Grade	G1	61	2.086	0.471, 1.697	0.101		6.070	1.953, 1.407	<0.001		4.153	1.758, 0.826	0.006
	G2	306		0.274, 2.280				2.706, 1.245				2.114, 1.243
	G3	125		0.696, 2.432				2.761, 1.456				2.351, 1.189
	G4	7		1.903, 1.194				2.801, 1.666				2.761, 0.602
Cancer stage	I	27	0.394	0.414, 1.766	0.853		0.714	2.407, 1.525	0.613		0.177	2.115, 1.178	0.971
	II	83		0.594, 2.271				2.669, 1.269				2.110, 1.403
	III	90		0.560, 2.146				2.401, 1.464				2.095, 1.129
	IVA	301		0.371, 2.359				2.645, 1.357				2.200, 1.177
	IVB	15		0.808, 1.920				2.800, 1.261				2.249, 1.143
	IVC	3		1.537, 1.242				3.028, 1.592				2.010, 1.889
Surgical margin status	Positive	59	0.631	0.714, 1.922	0.532		0.369	2.767, 1.178	0.692		−0.995	2.135, 1.229	0.371
	Close	50		0.316, 2.298				2.548, 1.184				2.302, 1.104
	Negative	349		0.374, 2.298				2.636, 1.418				2.052, 1.210
Perineural invasion	Yes	170	−2.145	0.069, 2.133	0.033		4.248	3.022, 1.303	<0.001		−5.032	1.762, 1.313	<0.001
	No	197		0.566, 2.306				2.431, 1.357				2.400, 1.081
Angiolymphatic invasion	Yes	122	0.412	0.411, 2.124	0.681		1.216	2.798, 1.408	0.225		0.620	2.169, 1.230	0.536
	No	232		0.309, 2.360				2.608, 1.374				2.084, 1.231

HNSCC, head and neck squamous cell carcinoma; HPV, human papillomavirus; SD, standard deviation.

Figure 1

Differential expression of miR-99a-3p, miR-411-5p and miR-4746-5p in HPV16+ HNSCC. (A) Volcano plot of the DEMs in HNSCC samples compared to normal tissues. (B) Volcano plot of DEMs in HPV16(+) HNSCC samples compared to HPV16(−) samples. (C) Venn diagram for common DEM selection. (D) Expression of miR-99a-3p, miR-411-5p, and miR-4746-5p in normal and HNSCC tissues. (E) Expression of miR-99a-3p, miR-411-5p and miR-4746-5p in HNSCC tissues with HPV16(+) and HPV16(−) status. (F) Expression of miR-99a-3p, miR-411-5p, and miR-4746-5p in HNSCC samples with different perineural invasion status. A violin plot displayed a median of log2 (miRNA-seq RPM expression +0.01) with 95% CI; unpaired t-test; *, P<0.05 was considered significant; ***, P<0.001. HNSCC, head and neck squamous cell carcinoma; DEMs, differentially expressed miRNAs; HPV, human papillomavirus.

Table 2

Univariate and multivariate Cox regression analyses in HNSCC patients

DEMs	Univariate analysis			Multivariate analysis
	HR (95% CI)	P value		HR (95% CI)	P value
miR-99a-3p	0.655 (0.481–0.891)	0.007		0.658 (0.483–0.898)	0.008
miR-411-5p	2.157 (1.337–3.482)	0.002		1.756 (1.018–3.029)	0.043
miR-432-5p	1.648 (1.206–2.252)	0.002		1.251 (0.846–1.849)	0.262
miR-495-3p	1.525 (1.121–2.075)	0.007		1.164 (0.807–1.677)	0.417
miR-758-5p	1.381 (1.013–1.885)	0.041		0.928 (0.646–1.333)	0.684
miR-4746-5p	0.680 (0.436–0.813)	0.001		0.650 (0.474–0.892)	0.008

HNSCC, head and neck squamous cell carcinoma; DEMs, differentially expressed miRNAs.

miR-99a-3p, miR-411-5p, and miR-4746-5p showed significant differential expression in HPV16+ HNSCC tumors

As shown in , both miR-99a-3p and miR-411-5p were frequently downregulated in HNSCC compared to normal samples (miR-99a-3p: 0.469±2.254 vs. 2.339±0.823, P<0.0001; miR-411-5p: 2.604±1.366 vs. 4.747±1.329, P<0.0001). The expression level of miR-4746-5p was significantly upregulated in HNSCC samples (2.159±1.204 vs. 0.368±1.019, P<0.0001). In HPV16+ HNSCC samples, miR-99a-3p and miR-4746-5p were significantly upregulated and miR-411-5p was downregulated [miR-99a-3p: HPV16(+) 1.263±2.435, HPV16(−) 0.306±1.977, P=0.037; miR-411-5p: HPV16(+) 1.614±1.429, HPV16(−) 2.739±1.320, P<0.001; miR-4746-5p: HPV16(+) 3.035±1.228, HPV16(−) 2.011±1.297, P<0.001, ].

Expression levels of miR-99a-3p, miR-411-5p, and miR-4746-5p reflecting PNI status of HNSCC patients

As shown in , the expression levels of miR-99a-3p, miR-411-5p, and miR-4746-5p in HNSCC patients with different clinical and pathological characteristics were also analyzed. The expression of all 3 miRNAs was significantly different with regard to PNI status. In HNSCC patients with positive PNI, miR-99a-3p and miR-4746-5p were significantly downregulated and miR-411-5p was upregulated (miR-99a-3p, P=0.033; miR-411-5p, P<0.001; miR-4746-5p, P<0.001, ).

Besides, miR-411-5p and miR-4746-5p showed significantly different expression levels in HNSCC patients with regard to cancer grade (miR-411-5p, P<0.001; miR-4746-5p, P=0.006), while miR-99a-3p showed significantly different expression in HNSCC patients with regard to cancer T stage (P=0.005). In addition, miR-99a-3p expression was upregulated in older HNSCC patients (≥60: 1.154±2.307, <60: 0.476±2.511, P=0.001). Furthermore, miR-411-5p expression was downregulated in HNSCC patients with a confirmed smoking history (smoker: 2.376±1.346, non-smoker: 2.676±1.368, P=0.034) and miR-4746-5p expression was upregulated in male HNSCC patients (male: 2.273±1.171, female: 1.867±1.508, P=0.001). All data are summarized in .

HNSCC patients with high miR-99a-3p and miR-4746-5p expression, along with low miR-411-5p expression displayed better OS

The OS of 518 HNSCC patients was analyzed in groups with different expression levels of miR-99a-3p, miR-411-5p and miR-4746-5p (). Cox regression analysis () revealed miR-411-5p as a risk factor (HR =1.756, P=0.043), and miR-99a-3p (HR=0.658, P=0.008) in addition to miR-4746-5p (HR =0.650, P=0.008) are protective factors in HNSCC. As shown in , HNSCC patients with low miR-411-5p (P=0.0013) and high miR-99a-3p expression (P=0.0068), along with miR-4746 (P<0.001) exhibited significantly improved OS. Furthermore, a risk model was constructed, based on the risk score of miR-99a-3p, miR-411-5p and miR-4746-5p, determined by multivariate Cox regression analysis. All HNSCC samples were categorized into a high-risk or low-risk group, based in the median risk score. HNSCC patients in the high-risk group showed significantly poorer OS than patients in the low-risk group (P<0.0001).

Figure 2

Overall survival of HNSCC patients in groups of negative/positive expression of miR-99a-3p (A), miR-411-5p (B) and miR-4746-5p (C) and high/low risk level (D). P<0.05 was considered significant. HNSCC, head and neck squamous cell carcinoma.

Figure 3

Functional enrichment analysis displaying significantly enriched GO processes (A) and KEGG pathways (B) of predicted target genes. GO, Gene Ontology.

Target gene prediction and functional analysis of miR-99a-3p, miR-411-5p, and miR-4746-5p

Three hundred and eighty-two predicted target genes of miR-99a-3p, miR-411-5p and miR-4746-5p were identified from the miRDB database. GO and KEGG pathway analysis were performed on Metascape with the cutoff P<0.05.

As shown in , significant BP terms included regulation of cellular catabolic process, gamete generation, regulation of protein catabolic process, developmental cell growth, neuroepithelial cell differentiation and neuron projection extension. The MF and CC terms were mainly enriched in signal sequence binding and the axon, respectively. Moreover, KEGG pathways were enriched in the MAPK and FoxO signaling pathways, which suggested that target genes may be associated with cancer and neural growth.

Correlation between the downregulation of epithelial mesenchymal transition (EMT) genes in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression

To determine the possible effects of miR-99a-3p, miR411-5p and miR-4746-5p in HPV16+ HNSCC patients, we screened the DEGs in HNSCC samples with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression prior to undertaking GSEA analysis. Analysis of hallmark gene sets displayed significant changes of gene expression related to EMT, angiogenesis, E2F-targets and myogenesis in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression (). EMT gene sets were significantly downregulated in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression ().

Figure 4

miR-99a-3p, miR-411-5p and miR-4746-5p deregulation were related EMT and invasion. (A) GSEA analysis for differently expression mRNAs in HNSCC tissues with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. (B) Significant down regulation of EMT gene set in HNSCC tissues with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. (C) Heat map and clustering of genes related to invasion in HNSCC tissues with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. (D) Comparison of genes related to invasion in HNSCC tissues with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. EMT, epithelial mesenchymal transition; GSEA, Gene Set Enrichment Analysis; HNSCC, head and neck squamous cell carcinoma.

Genes related to invasion were downregulated in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression

Selected genes related to invasion were compared, based on the significant relationship between 3 miRNA clusters with PNI status. As shown in , invasion-related genes (matrix metallopeptidase family members MMP1, MMP2, MMP3, MMP7, MMP9, MMP11, MMP12, MMP13, MMP16, and MMP28, integrin subunit beta 3 ITGB3, secreted protein acidic and cysteine rich SPARC) were significantly downregulated in HNSCC samples with miR-99a-3phighmiR-411-5plowmiR-4746-5phigh expression.

Discussion

High risk HPV are associated with precancerous lesions and carcinomas of cervical, vulvar, vaginal, anal, penile and HNSCC (18,19). HPV16 is the most common type in HNSCC and results in a distinct HNSCC subtype with unique pathological, clinical, molecular and prognostic characteristics (2,7,18,20).

Several miRNAs, including miR-21, miR-143, miR-9, miR-145, and miR-199a were identified in regulating the expression of host and viral factors that are particularly responsible for viral transcription which leads to HPV-related cervical carcinogenesis (21). HPV viral E6/E7 oncogenes may contribute to the deregulation of cellular miRNA network and subsequently affect tumor progression (22). miRNA deregulation, including miR-9, miR-21, miR-29a, miR-100 and miR-106b, were also reported in HPV16+ HNSCC (23-25). miR-15a/miR-16, miR-143/miR-145 and the miR-106/miR-363 cluster were suggested to target common HPV mediated pathways such as E6/p53 and E7/pRb in HPV-induced carcinogenesis (14). Furthermore, high risk E6 and E7 proteins could deregulate miRNA expression by directly interacting with other regulatory cellular proteins and DNA methylation (26,27). However, the interaction between HPV16 and the cellular miRNA network, along with the molecular mechanism behind tumor progression were not clear.

By comparing miRNA expression with the clinicopathological parameters of HNSCC patients, we revealed that less PNI events were identified in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh. PNI is well established as an independent prognostic indicator for HNSCC patients (28,29). In addition, HPV+ HNSCC patients were less likely to develop PNI than HPV- HNSCC patients (30). We proposed that miR-99a-3p, miR-411-5p and miR-4746-5p may participate in HPV16+ HNSCC progression and influence the patients’ prognosis.

To further determine the potential function of miR-99a-3p, miR-411-5p and miR-4746-5p, we explored the predicted target genes of 3 DEMs from the miRDB database and performed functional enrichment analysis. Based on the GO analysis, the predicted target genes of the 3 DEMs were enriched in cell growth, neuroepithelial cell differentiation, neuron projection extension, and neuron migration. KEGG pathway analysis showed that the predicted target genes were enriched in the MAPK signaling pathway and FoxO signaling pathway. The MAPK pathway is essential for the regulation of many cellular processes, including inflammation, cell proliferation, motility and apoptosis (31-33). It is well established that this pathway participates in HNSCC progression. Loss of FoxO function leads to increased cellular survival and plays a role in epithelial cancer progression (34).

We further screened DEGs in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. GSEA revealed genes related to the EMT pathway were significantly downregulated in HNSCC patients with miR-99a-3phigh miR-411-5plow miR-4746-5phigh expression. Additionally, genes related to invasion were also downregulated, including matrix metallopeptidase family members MMP1, MMP2, MMP3, MMP7, MMP9, MMP11, MMP12, MMP13, MMP16, and MMP28, integrin subunit beta 3 ITGB3, secreted protein acidic and cysteine rich SPARC (35-37). Overexpression of EMT related genes were widely reported in the perineural invasion (PNI) sites of many cancers, including pancreatic cancer, colorectal adenocarcinoma, salivary adenoid cystic carcinoma and HNSCC (38-41).

MMPs, expressed in cancer or stromal cells, are capable of degrading ECM to promote tumor invasion and metastasis (35,42). The overexpression of MMPs, including MMP1, MMP3, MMP9, MMP11, and MMP13, is considered a risk factor for PNI in several cancers (43-46). Saygili et al. showed that MMP2 activation, under frequency electrical field stimulation, contributed to pro-NGF conversion and induced nerve sprouting (47). The overexpression of MMPs, such as MMP2 and MMP9, has also been widely reported in HNSCC tissues. However, the molecular mechanisms in PNI have not been widely investigated (48-51). The present study indicated potential roles of miR-99a-3p, miR-411-5p and miR-4746-5p in EMT and PNI of HPV16+ HNSCC.

HPV+ HNSCC patients have displayed significantly more favorable outcomes compared to HPV- HNSCC patients? (52,53). In addition, HNSCC patients with different HPV statuses exhibit distinct clinical presentations, socioeconomic, and molecular profiles (6,54,55). However, predictive biomarkers for the outcomes specific to HPV-related HNSCC patients are not clear. miRNAs, which act as posttranscriptional gene regulators, have been shown to be involved in almost all key cellular processes in cancer progression. Tissue-specific miRNAs have been explored as potential biomarkers in many kinds of cancers (56-59). A variety of prognostic miRNA biomarkers have also been selected for HNSCC, such as let-7d, miR-205, miR-99b-3p, and miR-100a-5p (60,61).

In this study, DEMs were screened in HPV16-related HNSCC samples, based on the TCGA database. Through subsequent multivariate Cox regression analysis, miR-99a-3p, miR-411-5p and miR-4746-5p were selected as potential prognostic biomarkers for HPV16(−) HNSCC (). miR-99a-3p and miR-4746-5p were upregulated, and miR-411-5p was downregulated significantly in HNSCC tissues (). In addition, miR-99a-3p and miR-4746-5p predicted better outcomes, whereas miR-4746-5p forecasted poorer outcomes in HNSCC patients. As shown in , patients with miR-99a-3plow miR-411-5phigh miR-4746-5plow expression exhibited poorer OS. The dysregulation of 3 DEMs in HPV16+ cancer samples could serve as a predictor of improved prognosis in HPV16+ HNSCC patients. Previous investigations have demonstrated that miR-99a-3p is downregulated in prostate cancer, colorectal cancer and HNSCC (62-65). Furthermore, Okada et al. observed significantly suppressed cell migration and invasion of FaDu and SAS cells, post miR-99a-3p transfection (62). This demonstrated the potential antitumor function of miR-99a-3p and highlighted its potential as a predictive biomarker for the positive prognosis of HNSCC, consistent with the findings of this study. miR-411-5p is regarded as a tumor suppressor in non-small cell lung cancer, hepatocellular carcinoma, and breast cancer. However, its function in HNSCC has not yet been reported (66-68). Kumaran et al. explored breast tissue-specific expression of miR-4746-5p compared to adjacent normal tissue; however, no further study was conducted (69). In addition to miR-99a-3p, our study showed that miR-411-5p and miR-4746-5p could also serve as potential prognostic biomarkers for HPV16-related HNSCC.

Overall, we determined that miR-99a-3p, miR-411-5p and miR-4746-5p are expressed differently in HPV16+ HNSCC patients. HPV16+ patients showed less PNI and more favorable OS with miR-99a-3phighmiR-411-5plowmiR-4746-5phigh expression. Further functional analysis showed downregulation of EMT gene sets and invasion related genes in HNSCC patients with miR-99a-3phighmiR-411-5plowmiR-4746-5phigh expression. The above results indicated possible roles of miR-99a-3p, miR-411-5p and miR-4746-5p in HPV16+ HNSCC progression. Furthermore, more comprehensive studies are required to affirm the MFs of the 3 miRNAs in HPV+ HNSCC.

The article’s supplementary files as