Head and Neck Cancer Patients’ Survival According to HPV Status, miRNA Profiling, and Tumour Features—A Cohort Study

Head and neck cancers (HNC) are a heterogeneous group of tumours mainly associated with tobacco and alcohol use and human papillomavirus (HPV). Over 90% of all HNC are squamous cell carcinomas (HNSCC). Sample material from patients diagnosed with primary HNSCC (n = 76) treated with surgery as primary treatment at a single centre were assessed for HPV genotype, miR-9-5p, miR-21-3p, miR-29a-3p and miR-100-5p expression levels. Clinical and pathological data were collected from medical records. Patients were enrolled between 2015 and 2019 and followed-up until November 2022. Overall survival, disease-specific survival and disease-free survival were assessed and correlated with clinical, pathological, and molecular data. Kaplan-Meier and Cox proportional hazard regression was used to assess different risk factors. In the study, male gender, HPV-negative HNSCC (76.3%) mostly located in the oral region (78.9%) predominated. Most patients had stage IV cancer (47.4%), and the overall survival rate was 50%. HPV was found not to affect survival, indicating that in this population, classic risk factors predominate. The presence of both perineural and angioinvasion was strongly associated with survival in all analyses. Of all miRNAs assessed, only upregulation of miR-21 was consistently shown to be an independent predictor of poor prognosis and may thus serve as a prognostic biomarker in HNSCC.


Introduction
As a combined entity, head and neck cancer (HNC) is the eighth most common cancer globally with an incidence of more than 870,000 cases (lip, oral cavity, C00-C06; oropharynx, C09-C10; nasopharynx, C11; hypopharynx, C12-C13; larynx, C32) per year in the most recent GLOBOCAN 2020 data [1]. Approximately 90% of HNCs [2] are head and neck squamous cell carcinomas (HNSCCs), which derive from the mucosal epithelium in the upper aerodigestive tract [3]. Most HNSCCs, approximately 75%, share common risk factors including alcohol and tobacco use, while human papillomavirus (HPV) is the major risk factor for approximately 25% of HNSCCs [4]. Depending on the geographical region and culture, the mentioned risk factors have different importance for the incidence of HNSCC. For example, the decline in smoking prevalence in most high-income countries has facilitated the high prevalence of HNSCC in the US and Western Europe due to increased rates of oropharyngeal infection with carcinogenic HPV [3,5]. Thus, HPV-positive HNSCC is distinct from HPV-negative HNSCC regarding genetic, epigenetic, and protein expression profiles, epidemiological factors and clinical features, but the treatment for both groups is almost the same [6]. HPV is particularly associated with oropharyngeal squamous cell carcinoma (OPSCC), which affects the tonsils, base of the tongue, soft palate and uvula [5]. Also, HPV presence is an often-reported important factor when considering

Patient and Tumour Characteristics
Over the patient enrolment period (2015-2019), sample material was collected from 87 patients. However, in six cases, the material was not from the primary tumour; for three patients, the histopathologic diagnosis indicated that the tumour was not squamous cell carcinoma, and two patients had inoperable tumours. Thus 11 patients were excluded from further analyses.
According to clinical TNM (cTNM), most patients presented with stage IV cancer (47.4%). No patient had distant metastasis at the time of surgical treatment (M = 0); however, 17 patients developed metastasis during the follow-up (Table 1). Most patients with subsequent distant metastases presented with stage IV disease (11/17) and were HPV-negative (12/17). Clinically most patients were classified as T2, followed by T4a, T3 and T1. However, local lymph node involvement was seen in more than half of the patients (cN > 0, 53.9%, 41/76).
Samples were classified according to HPV status into two groups: 18 (23.7%) samples with HPV-positive and 58 samples with HPV-negative HNSCC (76.3%). The HPV-positive patients were only slightly younger than the HPV-negative group on average (59.6 ± 15.5 vs. 62.4 ± 10.1 years), and the difference was not statistically significant (t-test, p = 0.4749). Among HPV-positive tumours HPV16 was found in most cases (77.7%, 14/18). Almost 60% of patients underwent surgery as the only treatment, while the rest of the patients also subsequently received radiotherapy (Table 1). Only four patients received chemoradiotherapy, and they were grouped with those receiving radiotherapy and not reported separately.
Results of the histopathology assessment are presented in Table 2. Histopathological T stage followed a clinical pattern with T2 being the most represented, followed by T4a, T3 and T1. The concordance of cT and pT was high (Kappa interrater agreement 0.68; 95% CI 0.55-0.81), but it was lower for cN and pN classification (Kappa = 0.48; 95% CI 0.56-0.61). However, if patients with pNX were considered to be N0 since the tumour was considered small and surgical intervention minimal, then the concordance was again high (Kappa = 0.65; 95% CI 0.52-0.77). Pattern congruence also exists between clinical and pathological overall stage (Kappa = 0.74; 95% CI 0.61-0.87), with most staged identically while nine cases were up-staged on pathological classification and three down-staged. Histologically, most patients had grade II tumours (46.1%, 35/76). Angioinvasion or perineural invasion was present in the majority of tumours (65.8%, 50/76). There were no signs of any invasion in 26 patients. In the majority of patients (53.9%), more than 21 lymph nodes were excised. Half of the patients' lymph nodes were not involved ( Table 2). p-pathological, X-regional lymph nodes cannot be assessed, NS-not specified, NA-not applicable, LN-lymph node.

MiRNA Expression
Fold changes of assessed miRNA molecules (miR-9-5p, miR-21-3p, miR-29a-3p and miR-100-5p) compared to the normal sample pool across patient subgroups are shown in Table 3. Consistent upregulation of miR-21 can be seen across the study population and all subgroups ( Figure 1). MiR-9 was slightly more increased in HPV-positive subsets, while miR-29a-3p and miR-100-5p were downregulated across the groups. These findings are consistent with the HNSCC cohort of The Cancer Genome Atlas (TCGA) [18] (Supplementary Figure S1). Expression of other assessed miRNA molecules that could not be assessed in the complete sample population due to lack of material is presented in Supplementary Table S2. Table 3. Relative expression (median and IQR, fold change versus normal control sample) of miRNA molecules in the study population and subgroups.

Overall Survival
The median follow-up was 45.5 months, with a range from 0 to 82 months for the whole population. The median follow-up time of living patients was 61 months. Half of the patients were alive (38/76, 50%), and most of those living (36/38, 94.7%) had no evidence of disease, while 50% (38/76) of patients died ( Table 1). The overall survival rate was higher among women (70% vs. 42.9%). There is no considerable difference in survival rate between those with HPV-positive (10/18, 55.5%) and those with HPV-negative tumours (28/58, 48.3%). The overall survival rate in relation to the stage was as follows: 66.7% (14/21) stage I + II, 36.8% (7/19) stage III and 47.2% (17/36) stage IV.
The estimated five-year survival (Kaplan-Meier analysis) was 75.6% for stage I + II, 33.8% for stage III and 44.9% for stage IV. However, the difference between survival curves when staging was done according to clinical and pathological TNM variables was not statistically significant (p = 0.141 and p = 0.127, respectively) ( Figure 2).
Kaplan-Meier analysis was performed for all analysed parameters according to OS, DSS, and DFS survival and all curves are shown in Supplementary Figure S2.
Interestingly, only the presence of invasion ( Figure 3) and the high expression of miR-106b were consistently associated with significantly different survival curves for all three outcomes assessed (overall survival -OS, disease specific survival -DSS and disease free survival-DFS). In all cases finding both perineural and angioinvasion was associated with worse outcomes (OS p = 0.014; DSS p = 0.007; DFS p = 0.005; Figure 3). Not all patients could be assessed for miR-106b; however, for those available (n = 69), the survival curves were remarkably different (OS p = 0.020; DSS p = 0.024; DFS p = 0.026; Supplementary Figure S2). Notable is the lymph node yield analysis where the worst performing patients had 1-10 lymph nodes removed for analysis, while the best outcome was in those with 0 nodes dissected (OS p = 0.006; DSS p = 0.031; DFS p = 0.026). Clean resection margins were significant for DSS (p = 0.036) but not for OS (p = 0.138) (Supplementary Figure S2).

Risk Factors in HNSCC Patients
To assess the influence of clinical parameters as well as measured expression levels of miRNA molecules, univariate and multivariate Cox proportional hazard regression models for the OS, DSS and DFS were made. All univariate results are presented in Supplementary Table S3.
Multivariate models were created by entering variables with p ≤ 0.1 in univariate analysis, and the models are presented in Table 4. All three models were significant overall (OS, p < 0.001; DSS, p = 0.028; DFS, p < 0.001); however, for the DSS model, no variables were individually significant. For OS, gender, presence of perineural and angioinvasion, tumour resection margin, lymph node yield and the expression of miR-21 were independent prognostic factors. MiR-21 was also found to be significant in the DFS model.
Having noted the significant association of miR-21 in multivariate models as well as potential trends of miR-106b seen in the univariate analysis, we further investigated the effect of these miRNAs on survival across different stage categories using Kaplan-Meier analysis (Figure 4). High expression of either miRNA strongly affected the survival of patients with stage III or stage IV cancer, while the survival curves of stage I-II patients were almost not affected, especially for miR-106b (panels B, D-F green lines). Analysis of other miRNAs (-9, -29a and -100) did not yield significant differences and was not shown. * for the DSS model, only miR-106b was significant in univariate models, thus this miRNA was selected for multivariate models despite reducing the total number of cases evaluated to 69 due to missing data for some cases. HR-hazard ratio, CI-confidence interval, Sm-smoker, NSmND-never smoker/never drinker, SmD-smoker and drinker, LN-lymph node, FC-fold change, NA-not applicable. Association of individual miRNA expression (miR-9, miR-21, miR-29a and miR-100) with overall survival were evaluated with ROC curve analysis. However, the resulting AUC values were poor (between 0.52 and 0.57; Supplementary Figure S3).

Discussion
Over the last few decades, there has been substantial improvement in HNSCC understanding, diagnosis, staging and treatment [19]. The management of the patients is based on histologic parameters such as TNM staging and tumour grading. However, there is a need for new biomarkers in order to improve patients' prognosis and treatment [20].
HNSCC is one of the most aggressive cancers due to advanced disease at the time of diagnosis [13]. Most of the patients in the current study indeed presented with stage IV cancer (47.4%). Nevertheless, this number of stage IV patients' five-year mortality was 50%. The mortality rate of the patients was below the five-year European average agestandardized mortality reported in 2015 (65.6%) [21]. Although tobacco and alcohol were the predominant risk factors here, HPV is often also associated with HNSCC, especially in the case of OPSCC [22]. However, the majority of the tumours in this study were located in the oral region (78.9%), which is not the predilection site for HPV infection, and this might explain the diminished impact of HPV on survival in the present study. The HPV effect on OS in non-OPSCC is controversial [6], ranging from no prognostic benefit [23] to a good prognostic factor [24,25]. HPV16 was the most common type, which is in accordance with previous studies [5,26]; however, HPV-positive patients were only slightly younger in comparison with HPV-negative patients and the presence of HPV as determined by mRNA and DNA analysis was also not significantly associated with survival unlike in other more Western populations [27,28]. A similar non-existent HPV-associated survival benefit was observed in a previous retrospective analysis on an unrelated set of samples (n = 99 FFPE) [29], suggesting that different populations might have different contributions of HPV as a protective factor.
Approximately 40% of HNSCC were diagnosed in elderly patients (65+ years) ( Table 1), and this observation is in line with the current literature [21]. However, elderly patients frequently have comorbidities, which make them susceptible to other diseases and death [30], and this is also evident here, where more patients died from other causes (30.3%) than from HNSCC (19.7%) ( Table 1).
Our study confirmed earlier findings that the male population is more prone to both HPV-positive and HPV-negative HNSCC [31]. Possible reasons for such distribution of HPV-positive HNSCC could be associated with hormonal differences as well as different HPV transmission rates between genders, although the latter is still ambiguous according to the literature [31]. The observed gender inequality could also be explained by the fact that men are more likely than women to report increased numbers of sexual partners [5]. The culturally and traditionally conditioned presence of high-risk factors, such as tobacco and alcohol consumption, is more associated with the male population and could explain the predominance of males among those with HPV-negative HNSCC [3]. However, the uneven representation of women and men in this study could lead to potential bias. Thus, the weak significant association of gender with overall survival on multivariate Cox regression should be interpreted with caution (male vs female, HR 3.1 95% CI 1.1-8.7, p = 0.027, Table 4). There was no difference in stage at patient presentation between genders (chi square test, p = 0.962; Supplementary Table S1), thus, health awareness cannot be the cause for survival differences between gender. On the other hand, men were more subject to classic risk factors (smoking and drinking), and the difference was striking (15.0% female vs 75% men reporting alcohol and tobacco use, chi-square p < 0.001; Supplementary Table S1), which could easily explain the gender survival rate differences.
One significant prognostic factor was the presence or absence of perineural and angioinvasion. According to our observation, the worst survival could be observed in patients who had a combined type of invasion, while patients without any invasion had the best survival. The presence of invasion was strongly associated with survival in all versions of Kaplan-Meier and Cox analysis for all three examined outcomes (OS, DSS and DFS) suggesting a strong robustness of the association (Figure 3). Cox models also implied strong influence of lymph node yield where patients with 1-10 evaluated nodes had a relatively high risk of death in overall survival (HR 11.6 95% CI 1.2-9.7) and disease-free survival analysis (HR 10.0 95% CI 1.8-54.2) ( Table 4) compared with patients where no lymph nodes were dissected. This observation can be biased due to the extent of the disease, where small tumours are minimally excised and where lymph nodes are intentionally not evaluated.
Interestingly, clinical and pathological N stage or the actual number of positive lymph nodes found in the patient material or high LN ratio were never associated with survival on univariate Cox regression analysis (Supplementary Table S3) or Kaplan-Meier analysis (Supplementary Figure S2). Precise, numerical rather than a descriptive determination of the resection margins by the pathologist was implicated in survival (OS and DFS). Patients who received only a descriptive designation of resection margin ("clean") had a higher hazard ratio (OS HR 4.5, 95% CI 1.5-13.6) compared with patients where resection margin was 5 mm or more, which was even worse than the hazard ratio of patients with known close margins (OS HR 1.9, 95% CI 0.7-4.7) ( Table 4).
In accordance with some previous studies [13,15], our study showed that elevated miR-21 levels were an independent predictor of poor prognosis (Table 4), with especially strong differences in stage IV cancer ( Figure 4A-C). Depending on the tumour type, anatomic site or cellular context, miR-9 can act as an oncomiR or oncosuppressor [13]. In this study, we observed an increasing trend of miR-9-5p in HPV-positive tumours (Figure 1). However, possibly due to the limited number of samples from different sub-locations, any effects on survival were not significant. Additionally, a high expression level of miR-100 was correlated with poorer overall survival (Supplementary Figure S2), similar to some other studies [18,32]. However, miR-100 did not reach statistical significance in Kaplan-Meier or Cox proportional regression analysis. Although it was not possible to detect miR-106b in the whole set of samples, decreased miR-106b was correlated to better survival, and this finding was statistically significant in some cases of univariate Cox analysis (Supplementary Table S3) as well as in Kaplan-Meier analysis when combined with tumour stage (Figure 4D-F).
The limitations of our study include the relatively small number of patients and potential but unavoidable gender bias. Other limitations include the possible reporting biases and underreporting of important cancer features as evidenced in the variability of pathologists' reports, with some describing resection margins as narrow or wide rather than providing precise measurements in some cases, some explicitly noting different invasions or the presence of ENE. Due to the inconsistencies in reporting the depth of invasion, it was completely impossible to assess this as a separate parameter in most cases. Unfortunately, almost no p16 staining was performed, except in 4 out of 76 cases, limiting the assessment of this useful biomarker in our study population. The implementation and completion of a uniform report for HNC would reduce this issue.

Patients and Tumour Samples
All patients in this cohort study underwent surgery as primary treatment for HNSCC performed at the Department of Maxillofacial Surgery, Clinical Hospital Dubrava, University of Zagreb, Croatia, between 2015 and 2019 and followed-up until November 2022. Patients with a diagnosis of oral cancer (O; gingiva, retromolar area, oral tongue, sublingual area-excluding base of the tongue, buccal mucosa), as well as oropharyngeal cancer (OP; base of the tongue, tonsil, posterior pharyngeal wall) were included. In our previous study [33], a subset of patients was assessed by high-throughput methods for miRNA profiling. However, at that time, the patient follow-up was too short to make observations regarding outcomes. Within the current study, we collected follow-up information and included additional patients treated in the intervening period.
A total of 76 patients, 20 women (age range 32-87 years) and 56 men (age range 31-85 years), were included in the current study. Patients' data from previously collected (n = 59) and newly enrolled patients (n = 17) are shown in Supplementary Table S1. For the current study additional detailed information including clinical tumour, node, and metastasis (cTNM) status, pathological TNM (pTNM) status, histological grade, presence of histopathologically assessed angioinvasion or perineural invasion or their combination, tumour involvement of surgical margins (and the distance to the margin), lymph node yield (LNY), positivity and ratio (LNR), presence of extranodal extension (ENE), postoperative treatments and survival information including disease recurrence (Tables 1 and 2) were collected from hospital databases and medical records for all included patients. Tumours were staged by using the 8th edition of the AJCC Cancer Staging Manual [9]. Since p16 information was not available, cases were staged according to the p16 negative guidelines. Separate staging is shown for clinical and pathological TNM classification (Tables 1 and 2). Furthermore, alcohol and tobacco consumption as risk factors were reviewed; however, pack/year or detailed alcohol consumption data were not available.
Informed consent was obtained from all patients, and the study (Epic-HNSCC project No 4758) was approved by the Clinical Hospital Dubrava Bioethics Committee (EP-KBD-10.06.2014) and the Ruder Bošković Institute Bioethics Committee (BEP-3748/2-2014).

RNA and DNA Extraction and HPV Testing
Detailed procedures were previously described [33] and were also applied to the newly enrolled patients. Briefly, DNA and RNA were isolated from two separate tumour tissue samples using an EZ1 DNA Tissue kit (Qiagen, Hilden, Germany), GenElute-E Single Spin Tissue DNA Kit (Sigma-Aldrich, St. Louis, MO, USA) and miRNeasy Mini Kit (Qiagen, Hilden, Germany) for the isolation of DNA and RNA, respectively. The total RNA and DNA concentration of samples was measured using a NanoPhotometer (Implen, München, Germany).
Polymerase chain reaction (PCR) for HPV was performed with consensus HPV-specific primers PGMY, GP5+/GP6+, LC, HPV16 and HPV18 type-specific primers as described previously [33] and detailed in the Supplementary methods document. Beta-globin PCR amplification was used as an internal control. The E6 mRNA analysis was performed on HPV16 DNA-positive samples as previously described [33]. Only HPV16 mRNA-positive and HPV18-positive samples were considered likely to be HPV associated.

MicroRNA Quantitation
The expression of miRNA targets was assessed using TaqMan advanced miRNA Assays (Applied Biosystems, Waltham, MA, USA) as described before [33]. cDNA was prepared from total RNA from samples using the TaqMan Advanced miRNA cDNA Synthesis Kit (Applied Biosystems, Waltham, MA, USA). Expression of miR-9-5p, miR-21-3p, miR-29a-3p and miR-100-5p were selected for further investigation in the current study. miR expression was normalised using the average of miR-16-5p, miR-181-5p and miR-191-5p as an internal reference control. Quantitative RT-PCR reaction was performed using CFX96 Touch Real-Time PCR Detection System (BioRad, Hercules, CA, USA) for newly collected samples. Cycle conditions and reaction volumes used for poly(A) tailing, adaptor ligation, reverse transcription, miR-Amp reaction and qPCR were according to the manufacturer's instructions. A pool of cDNA from three healthy tonsil samples was used as the referent sample [33]. The 2 −∆∆Ct method was used to calculate the fold changes. Due to a lack of material in some cases miR-106b-5p, miR-143-3p, miR-145-5p, and miR-199b could not be analysed for all samples, but old data were re-evaluated in the context of additional collected clinical, as well as survival data. Furthermore, to avoid issues with proper controls for normalization of the 2 −∆∆Ct calculations, obtained expression data were also analysed as single 2 −∆Ct versus internal reference miRNAs only. This additional 2 −∆Ct relative expression analysis is only presented as supplementary material.

Statistical Analysis
All statistical analyses were performed using the MedCalc software (version 20.111). Patients were dichotomised by the median expression of each studied miRNA into above and below median groups. Groups were also divided into above or below 65 years of age. Variables with levels consisting of few cases were grouped, i.e., stage I (n = 3) and stage II (n = 18) cancers were combined, as well as different sublevels of stage IV (IVa, n = 32 and IVb, n = 4) to allow more meaningful comparisons. Overall survival (OS-patients who died of any cause), disease-specific survival (DSS-patients who died of the main disease) and disease-free survival (DFS-any recurrence, distant metastasis or death of any reason) were recorded. Follow-up time in months was calculated from the date of enrolment (surgery) to the registered time of death or last registered follow-up in case of OS or DSS. Follow-up time in months from enrolment to disease recurrence, metastasis or death of any cause was considered for disease-free survival assessment. Individual parameters were assessed for influence on survival using the Kaplan-Meier method and log-rank test. Cox proportional hazard regression was used to assess the impact of examined parameters on the outcomes using the same time and event criteria as for the Kaplan-Meier method. Both univariate and multivariate models, with relevant variables, were constructed. Final multivariate models were constructed using non-redundant variables with p < 0.1 in univariate analysis (i.e., if both clinical and pathological TNM stage variables had p < 0.1, only the clinical one was included in the multivariate model). The performance of different miRNAs was also analysed with ROC curve analysis for overall survival. All p-values below 0.05 were considered statistically significant.

Conclusions
HNSCC is a highly aggressive cancer with nearly 50% of patients having the advanced disease at the time of diagnosis and a 50% mortality rate. Tobacco and alcohol continue to play a leading role as risk factors in HNSCC within the Croatian population, especially among men, while HPV seems to have a lower impact. However, a wider sample panel and p16 testing would be useful to draw a definitive conclusion. On the other hand, our results indicate that certain miRNAs, such as miR-21 and miR-106b, might serve as useful prognostic biomarkers in HNSCC, in addition to tumour features like perineural infiltration, angioinvasion and lymph node yield.

Data Availability Statement:
The data presented in this study are available in the Supplementary  Table S1.

Conflicts of Interest:
The authors declare no conflict of interest.