miR-154 Influences HNSCC Development and Progression through Regulation of the Epithelial-to-Mesenchymal Transition Process and Could Be Used as a Potential Biomarker

MicroRNAs and their role in cancer have been extensively studied for the past decade. Here, we analyzed the biological role and diagnostic potential of miR-154-5p and miR-154-3p in head and neck squamous cell carcinoma (HNSCC). miRNA expression analyses were performed using The Cancer Genome Atlas (TCGA) data accessed from cBioPortal, UALCAN, Santa Cruz University, and Gene Expression Omnibus (GEO). The expression data were correlated with clinicopathological parameters. The functional enrichment was assessed with Gene Set Enrichment Analysis (GSEA). The immunological profiles were assessed using the ESTIMATE tool and RNAseq data from TCGA. All statistical analyses were performed with GraphPad Prism and Statistica. The study showed that both miR-154-5p and miR-154-3p were downregulated in the HNSCC samples and their expression levels correlated with tumor localization, overall survival, cancer stage, tumor grade, and HPV p16 status. GSEA indicated that individuals with the increased levels of miR-154 had upregulated AKT-MTOR, CYCLIN D1, KRAS, EIF4E, RB, ATM, and EMT gene sets. Finally, the elevated miR-154 expression correlated with better immune response. This study showed that miR-154 is highly involved in HNSCC pathogenesis, invasion, and immune response. The implementation of miR-154 as a biomarker may improve the effectiveness of HNSCC treatment.


Introduction
Head and neck squamous cell carcinomas (HNSCCs) account for about 90% of all cancers in the head and neck area and are among the most common neoplasms worldwide [1]. This group includes tumors of the upper part of the digestive and respiratory system, particularly the oral cavity, larynx, and pharynx [2]. HNSCC is typically diagnosed at the advanced stage, resulting in high morbidity and mortality [3]. It is mainly associated with exposure to tobacco, alcohol abuse, and human papillomavirus (HPV), particularly type 16, which plays a crucial role in HNSCC development [4,5]. HPV-positive HNSCC tumors are more frequently found in the oropharyngeal area, with a younger age of onset, and also show less molecular diversity and higher sensitivity to radiation or chemotherapy when compared to HPV-negative tumors [6]. Moreover, HPV-positive HNSCC patients have significantly longer overall survival (OS) [6]. However, the five-year survival rate in HNSCC is significantly lower compared with other malignancies such as breast, cervical, or colorectal cancer [7]. The high mortality rate is associated with a lack of diagnostic biomarkers, leading to failure in early diagnosis and insufficient effectiveness of therapeutic methods [7]. Molecular biomarkers are important tools to diagnose, assess the likely course of the disease, and predict response to the treatment. Additionally, they are necessary to implement personalized treatment [8,9]. Although many biomarkers seem to affect the diagnosis and prognosis of HNSCC patients, only a few of them are approved for clinical use [7,10].
MicroRNAs (miRNAs) are single-stranded ncRNA molecules with 21-23 nucleotides in length. They are involved in the post-transcriptional regulation of gene expression through translational inhibition or even destabilization of the targeted messenger RNA (mRNA) [11]. They play a pivotal role in maintaining tissue homeostasis by regulating processes such as cell proliferation, differentiation, or apoptosis [12]. miRNAs can act as tumor suppressors or oncogenes. The discovery of their exact role is pivotal for understanding the molecular pathways involved in the development of HNSCC, both HPV-negative and HPV-positive [13][14][15][16]. miRNAs have been extensively studied as potential cancer biomarkers mainly due to their stability in a variety of biological materials, their specificity, and because they are easy to process and analyze with current protocols [17].
In this study, using data from The Cancer Genome Atlas (TCGA), we analyzed the expression and biological role of miR-154 strands 5p and 3p and determined their potential utility as HNSCC biomarkers. It is known that the expression level of miR-154 is downregulated in melanoma [18], glioma [19], non-small-cell lung (NSCLC) [20], breast [21], bladder [22], gastric cancers [23], and laryngeal squamous cell carcinoma [24] and upregulated in renal cell carcinoma [25]. Analysis of miR-154 functions revealed that this miRNA acts as a regulator of proliferation, cell viability, invasion, and migration, apoptosis of cancer cells in vitro and regulates tumor growth in vivo [18,[22][23][24][25][26]. Nevertheless, there is a lack of studies discussing how each miR-154 strand is expressed and how they affect the molecular characteristics of the tissue or the immune profile of the patients. In recent years, many reports have emphasized the importance of analyzing both miRNA strands separately [27][28][29][30]. It has been proven that these molecular species can be co-expressed in different tissues and display different modulatory roles in tumorigenic processes independently and/or cooperatively [27][28][29][30]. Therefore, within this study, we performed such an analysis of miR-154 for the first time. Our in silico investigation was based on a cohort of HNSCC patients and aimed to elucidate whether miR-154 could be further considered for clinical applications.
Additional information about the tumors' biological and genetic features was obtained from the supporting data presented by Thorsson et al. [36].
All the data are available online with unrestricted access and does not require the patients' consent or other permissions. The use of the data does not violate the rights of any person or any institution.

Functional Enrichment Analysis and Prediction of Genes' Function
Gene Set Enrichment Analysis (GSEA) software version 3.0 (http://www.gsea-msigdb. org/gsea/index.jsp (University of California, San Diego, USA, and Broad Institute, USA, accessed on 9 December 2021) was used for the analysis of functional enrichment [37]. The HNSCC patients were divided into two groups with high and low expression of miR-154 using the average expression as a cutoff. The input file contained expression data for 20,530 genes and 512 patients. We used 1000 gene set permutations for the analysis and pathways (oncogenic signatures (C), hallmark gene sets (H), and collection from MSigDB) with nominal p-value ≤0.05 and FDR q-value ≤0.25 considered significant.

Immune Profile Analysis
Analysis of the ESTIMATE, immune and stromal scores (Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data) required downloading a specific dataset from https://bioinformatics.mdanderson.org/estimate/disease.html (MD Anderson Cancer Center, USA, accessed on 9 December 2021) [47]. Then, those scores were used to define the infiltration of immune cells into tumor tissues and infer tumor purity. Subpopulations of specific immune cells were assessed using supporting data presented by Thorsson et al. and compared between subgroups of the patients with high and low expression of miR-154-5p, as well as of miR-154-3p [36]. Individuals were assigned into groups based on the mean value of the studied miRNA expression level.
Populations of specific immune cells, lymphocytes (T and B cells with their subpopulations), neutrophils, eosinophils, mast cells, dendritic cells, and macrophages (M0, M1, and M2), as well as T cell receptor (TCR) and B cell receptor (BCR) parameters, TCGA immune subtypes (C1-C4 and C6), along with other immune features such as IFN gamma and transforming growth factor beta (TGF beta) response, macrophage regulation, lymphocyte infiltration, TIL regional fraction, intratumor heterogeneity, proliferation, wound healing, homologous recombination defects, SNV and indel neoantigens, silent and non-silent mutation rates, aneuploidy, and cancer/testis antigen score (CTA score) were determined by conducting analyses of the supporting data presented by Thorsson et al. [36]. The data downloaded for this study can be explored and visualized through the interactive portal CRI iAtlas (www.cri-iatlas.org (Institute for Systems Biology, USA, accessed on 9 December 2021)) [48]. The analysis was performed as described in Section 2.6 "Statistical Analysis".

Statistical Analysis
All the statistical analyses were performed using GraphPad Prism 8 (GraphPad, San Diego, CA, USA) and Statistica v13.1 (Dell Software, Round Rock, TX, USA). The Shapiro-Wilk normality test, t-test, or Mann-Whitney U test was used for miR-154-5p and miR-154-3p expression levels (depending on the clinical parameters). The level of miR-154-5p and miR-154-3p (depending on the National Institute of Health classification of cancer location, as well as the TCGA molecular subtype) was determined using one-way ANOVA with Tukey's multiple comparisons test or the Kruskal-Wallis test with Dunn's multiple comparisons test. A negative correlation between miR-154-5p and miR-154-3p expression levels, as well as protein-coding genes, were analyzed in Statistica (correlation matrix) or GraphPad Prism 8 with the Pearson or Spearman's correlation tests depending on the data distribution. All the TCGA data are displayed as the means with the standard error of the mean (SEM). For DSF and OS analyses, the log-rank (Mantel-Cox) and Gehan-Breslow-Wilcoxon tests were used, and the hazard ratio (Mantel-Haenszel; HR) and the 95% confidence interval (CI) of the ratio were calculated. Graphs presenting boxes use 5-95 percentile whiskers, and column bars show the means with 95% CI. In all the analyses, p < 0.05 was used to determine statistical significance.

Availability of Data and Materials
The datasets used and/or analyzed during this study are available from the corresponding author on reasonable request. Raw data are available on the cBioPortal, in the UALCAN and the University of California, Santa Cruz databases, and the published Supplementary Materials.

Results
miR-154-5p and miR-154-3p were downregulated in the cancer samples when compared to the cancer-free samples, in pharyngeal squamous cell carcinomas and in the atypical molecular subtype of HNSCC.

miR-154-5p and miR-154-3p Levels Differ Depending on the Clinicopathological Parameters
The expression levels of miR-154-5p and miR-154-3p were analyzed in the previously determined subgroups depending on the available clinicopathological parameters for all the HNSCC samples. Significant differences in the expression levels for miR-154-5p were observed in the N0 vs. N1 + N2 + N3 cancer N-stages (2.071 ± 0.8426 vs. 1.866 ± 0.8409, p = 0.0066), I + II vs. III + IV cancer stages (2.145 ± 0.8002 vs. 1.913 ± 0.8548, p = 0.0077), G1 + G2 vs. G3 + G4 cancer grades (1.913 ± 0.8247 vs. 2.182 ± 0.8299, p = 0.0017), and in the negative HPV p16/ish samples in comparison to the positive ones (1.492 ± 0.9784 vs. 1.943 ± 0.7552, p = 0.0102). The other analyzed parameters did not differ between the studied groups for miR-154-5p. Significant differences between the miR-154-3p expression levels were detected only between the N0 vs. N1 + N2 + N3 cancer N-stages (1.623 ± 0.8086 vs. 1.455 ± 0.7938, p = 0.0208) and in the negative HPV p16/ish samples in comparison to the positive ones (1.208 ± 0.7566 vs. 1.534 ± 0.6708, p = 0.0277); the other analyzed parameters did not vary between the studied groups. All the data are presented in Table 1. Next, the results based on the TCGA data were tested using datasets GSE31277 and GSE144711 from the GEO database, and the expression levels of miR-154-5p and miR-154-3p depending on sample type, stages of cancer, and histological differentiation were assessed. In the case of GSE31277, significantly lower expression levels of miR-154-5p in the cancer samples in comparison to the surgical margin were observed (10.58 ± 0.1930 vs. 11.39 ± 0.1546, p = 0.0028), and no differences were observed for miR-154-3p (10.79 ± 0.2095 vs. 10.54 ± 0.1884, p = 0.3840). Moreover, no differences in the expression levels between the T-and N-stages as well as histological differentiation for both miR-154-5p ( Figure 2B).

Patients with Low Expression of miR-154-3p Had a Significantly Extended Overall Survival
The HNSCC patients were divided into low and high miR-154-3p and miR-154-5p expression groups using the mean expression as a cutoff, and then the 5-year OS and DFS were assessed. Slight changes in the OS time between low and high miR-154-5p expression groups were observed (p = 0.0706 and p = 0.0905, respectively), and no differences (p > 0.05) for the DFS were seen ( Figure 3A). However, a significantly longer 5-year OS was discovered in the patients with low miR-154-3p expression levels in comparison to the high expression group (p = 0.0013 and p = 0.0094, respectively). In the case of DFS, no difference between the patients with low and high expression levels of miR-154-3p was determined ( Figure 3B).

Higher Expression of miR-154-5p and miR-154-3p Was Associated with the Upregulation of Oncogenic Pathways in the HNSCC Patients
The biological influence of high and low expression levels of miR-154-5p and miR-154-3p was analyzed using GESA. In the case of all the HNSCC subsites, significant enrichment scores (p < 0.05 and FDR q-value <0.25) were observed for 36 gene sets from oncogenic signatures and five gene sets from cancer hallmark signatures for the patients with higher levels of miR-154-5p compared to the ones with lower expression of this miRNA. Similarly, for high levels of miR-154-3p, significant enrichment in 41 gene sets from oncogenic signatures and two gene sets from cancer hallmark signatures were detected. In both cases, those genes were implicated in AKT-MTOR, CYCLIN D1, KRAS, EIF4E, RB, as well as ATM pathways. Moreover, patients with higher levels of both strands displayed gene expression signatures characteristic of cancer stem-like cells, with features including genes downregulated during early stages of differentiation or connected with the WNT pathway. The above was supported by the enrichment of 194 genes associated with the Epithelial-to-Mesenchymal (EMT) process (NES = 1.930 for miR-154-5p and NES = 1.813 for miR-154-3p, respectively). It should be noted that gene enrichment analysis depending on the cancer subsite (oral cavity, larynx, or pharynx) also suggests a strong association of higher levels of both miR-154 strands with an aggressive phenotype of HNSCC. All the data are presented in Figure 4A,B and in Supplementary Table S1. Next, gene sets from the processes selected by means of GSEA were screened for potential targets for miR-154-5p using miRNA target prediction databases, and the expression correlation between miR-154-5p and the genes was calculated. Based on this approach, 24 genes were selected: VPS4B, MYO6, LRRC1, DMXL1, CEACAM7, OVOL1, ABI1, CDKN2B, NPEPPS, ESRP1, CPEB3, XPO7, LIN54, CNOT4, PLAGL2, C12ORF29, PTBP3, TMPRSS11D, WAC, LIN9, RNF138, GNAI3, DCAF16, and NSD2 (with R-coefficient from −0.09532 to −0.3055 and p < 0.05) ( Figure 5A). Using the GeneMANIA tool, the network of interactions and functions between the miR-154-5p targets was indicated. It was observed that those genes were co-expressed (61.19% of genes), had physical interactions (24.85% of genes), shared the same colocalization (7.83% of genes), and were from the same pathway (5.28% of genes). Moreover, these 24 genes are involved in the endosomal sorting complex required for transport (ESCRT), cell cycle, and viral replication ( Figure 5B). Finally, the association of the miR-154-5p targets with the HNSCC patients' survival was determined. In the case of higher expression of CEACAM7 and PLAGL2, a significantly longer DFS was observed (HR = 0.58, p = 0.022; HR = p = 0.036). Surprisingly, only for higher expression levels of CPEB3, longer OS time was observed (HR = 0.56, p = 0.0045). However, when the expression levels of all the 24 genes were compiled, a strong association between higher expression levels and longer OS time for the HNSCC patients was observed (HR = 0.58, p = 0.0051). No differences for DFS time were observed (HR = 0.71, p = 0.13) ( Figure 5C). The detailed data are presented in Supplementary Table S2.

Patients with Low Levels of Both miR-154 Mature Strands Have a More Favorable Immunological Profile
The last step of this study was to assess the immune profile of the patients with high and low expression levels of miR-154-5p and miR-154-3p. The ESTIMATE analysis of the previously determined patient groups distinguished based on the mean miRNA expression values was conducted, and then the lymphocyte infiltration score and the distribution of the immune cell fractions among them were examined. Investigation of the ESTIMATE score, which determines the cellular composition (purity) of the tumor, showed statistically significant differences between individuals with high and low levels of miR-154-5p (p = 0.0112) as well as of miR-154-3p (p = 0.0373). Next, the immune score was analyzed, but no significant (p < 0.05) differences between the studied subgroups were observed. Subsequently, the stromal score, which indicates the presence of stroma in tumor tissue was assessed. It was indicated that the patients with low expression of both miR-154 strands had a significantly decreased stromal score (p < 0.0001 and p < 0.0001, respectively). All the data are presented in Figure 6A.
Even though the immune score did not show any differences between the studied groups, the infiltration of lymphocytes varied. Individuals with low expression levels of both miR-154-5p and miR-154-3p showed a higher lymphocyte infiltration signature score compared to the high expression subgroups (p = 0.0377 and p = 0.0014, respectively) ( Figure 6B).

Figure 7.
Immunological characteristics including the estimated CTA score, TGF beta response, SNV neoantigens, and wound healing (A) and the level of silent and non-silent mutations, as well as homologous recombination defects (B) depending on high and low expression of miR-154-5p and miR-154-3p in the HNSCC patients; t-test, p < 0.05 considered as significant; ns-not significant, * p < 0.05, ** p ≤ 0.01, *** p ≤ 0.001.

Discussion
Head and neck cancers are a group of neoplasms with a strong ability to metastasize to local lymph nodes, which in consequence leads to a high mortality rate [1,49]. Moreover, the standard treatment based on radiotherapy, chemotherapy, as well as immunotherapy does not bring satisfactory results and needs to be improved [50][51][52]. Personalization of therapy requires the use of biomarkers, which could help to diagnose and predict the patient outcome. The molecular biomarkers seem to be the best option because they may describe molecular states of the cell and dictate treatment strategies [53][54][55]. Among the widely studied types of molecules, which could become valuable oncology biomarkers, are RNAs. It was demonstrated that both protein-coding and non-coding RNAs such as miRNAs or long non-coding RNAs (lncRNAs) are dysregulated in HNSCC and are potentially useful in clinical practice [56][57][58][59][60][61].
In this study, we used data from TCGA to analyze the expression and the biological role of miR-154-5p and -3p strands, as well as their usefulness as biomarkers in HNSCCs.
Our first important observation was that miR-154-5p and miR-154-3p are downregulated in HNSCC patients, and miR-154-5p specifically distinguishes cancer-free from cancer tissue. GEO datasets were used to validate the results obtained with TCGA data. Our results corroborate the previous studies regarding laryngeal squamous cell carcinoma, showing that the expression level of miR-154 in tumors is decreased in comparison to the healthy tissue [18][19][20][21][22][23][24].
It was observed that the expression levels of miR-154 are significantly downregulated in oral squamous cell carcinoma (OSCC) tissues in comparison to tumor-free surgical margins [34].
Moreover, we indicated that neoplasms situated in the pharynx had the lowest expression levels of miR-154; however, the difference was statistically significant only when compared to the oral cavity. Niu et al. analyzed the expression of miR-154 and demon-strated its downregulation in laryngeal cancers, although the study focused only on the Chinese population [24]. In fact, extensive research concerning different HNSCC subsites is lacking. Our study analyzed TCGA data obtained from samples collected predominantly from Caucasians, and not Hispanic or Latino populations, but the results were confirmed using a GEO dataset generated from individuals from South America, improving the diversity and relevance of the findings.
We found that low expression of miR-154-5p and miR-154-3p was significantly associated with a higher N-stage. Decreased miR-154-5p was also linked with the III + IV cancer stages. Surprisingly, our results obtained using the GEO database were opposite; the expression of miR-154-5p was higher in stages III + IV. In agreement with TCGA dataset results, other studies showed that low levels of miR-154 correlate with more advanced disease (higher grade and/or cancer stages) [24,62,63] and could reflect different study populations. We observed the trend of higher expression levels of miR-154-5p depending on the histological differentiation status, where expression of that miRNA was the highest in the poorly differentiated cancers (GEO) and within the mesenchymal subtype (TCGA). However, as it had been previously indicated in gliomas, low expression of miR-154 was linked to the higher grade, large tumor size, and unfavorable Karnofsky score [19]. The same situation was described in the patients with bladder cancer, where miR-154 downregulation was significantly associated with more advanced tumor staging, higher histologic grades, and lymph node metastasis [22]. On the other hand, in melanoma tissues, the increased level of miR-154 correlated with higher stages of tumor development and ulceration [18].
In this study, there were no significant associations between miR-154 expression and age or sex in the patients from our dataset, which stays in agreement with other studies on various cancers [18,19,22]. Furthermore, we observed that the patients with positive HPV p16 status displayed low levels of both miR-154 strands. Vojtechova et al. indicated that two members of the miR-154 family, miR-323a-3p and miR-487a-3p, were dysregulated in HPV-positive compared to HPV-negative tonsillar tumors [64]. miR-154-5p and -3p are downregulated in HPV-positive cervical cancer cell lines, but their roles were not examined by the authors [65]. This observation was confirmed by Zhao et al. using different cervical tissues and cell lines. Knockdown of miR-154-5p in HPV-16-positive SiHa cells caused increased cell proliferation, migration, along with the invasive ability, and its upregulation led to the reversed effect. They proved that miR-154-5p directly binds to the 3 UTR of CUL2, the core component of the E3 ubiquitin protein ligase complex, which, in turn, regulates the steady-state levels and stability of pRb [66]. It is worth mentioning that no reports of the direct role of miR-154-5p or -3p in HPV-positive HNSCC are available, and their diagnostic and biological role should be determined in the future. The ambiguity of associations between miR-154 strands and various clinicopathological parameters corroborates the complex nature of this regulatory system. To draw precise conclusions, complex interaction networks should be considered, which we attempted with GSEA.
We investigated whether miR-154-5p and miR-154-3p could be used as prognostic biomarkers for the assessment of individuals with HNSCC. We observed that the patients with lower levels of miR-154-3p and miR-154-5p had better OS compared to the group with high expression. The results from a study by Lin et al. confirmed this finding in renal cell carcinoma; however, the expression level of this miRNA was upregulated in cancer tissues [25]. Niu et al. described the reverse correlation between miR-154 and laryngeal cancer patient survival compared to our results [24]. Moreover, the same observation was proved for melanoma [18], glioma [19], and bladder cancer [22]. However, there is no available information about the influence of the miR-154 expression level on the OS in NSCLC [20], breast cancer [21], and gastric cancer [23]. It should be emphasized that there are no other independent studies concerning miR-154 in HNSCC than that presented by Niu et al. The authors were aware of the discrepancies between the results of their analyses, e.g., of the association of low miR-154-3p with a higher N-stage (tumorsuppressive character) and better OS (oncogenic role); however, these types of correlations rarely provide unequivocal results when based on a single miRNA. Moreover, it is well-known that the biology of neoplastic cells differs significantly from that of normal cells, and primary neutral miRNAs regulating tens or even hundreds of mRNA transcripts can be turned into pro-or antitumor molecules [67].
GSEA showed an association of higher levels of miR-154-5p and miR-154-3p with the enrichment of genes from oncogenic and cancer hallmark signature sets. The most dysregulated genes were related with the signaling pathways such as AKT-MTOR, CYCLIN D1, KRAS, EIF4E, RB, and ATM, which play an essential role in cancer cell proliferation, survival, and response to external stimuli. For example, E2F transcriptional activity is tightly regulated throughout the cell cycle via transcriptional and translational regulation, post-translational modifications, protein degradation, binding to cofactors, and subcellular localization. Furthermore, E2F participates in apoptosis and cell differentiation. Alterations in components of this pathway coincide with poor prognosis in cancers, emphasizing their importance for the clinical cancer phenotype. An intriguing addition to the understanding of E2F crosstalks was the finding that their activity could be regulated by miRNAs whose dysregulation has been implicated in malignancy. In turn, miRNAs themselves are targets of the E2F family proteins establishing negative feedback loops [68,69]. In breast cancer tissues, E2F5 was identified as a target of miR-154 with the inversely correlated expression. Additionally, it has been reported that in breast cancer cell lines (MCF-7), overexpression of miR-154 significantly inhibited cell proliferation, migration, and invasion, as well as increased cell arrest at the G0/G1 stage. These findings indicate that miR-154 acts as a tumor suppressor targeting E2F5 in breast cancer [70]. In this study, we proved that genes from the E2F target dataset had a decreased expression in the group of patients with high miR-154 levels. Moreover, those genes are potential targets of miR-154-5p. Thus, E2F could be recognized as a potential target of miR-154, and as a result of its degradation, there may be a reduction in the expression of these transcription factors' target molecules.
The genes involved in the G2/M checkpoint play a crucial role in progression through the cell division cycle. The G2 checkpoint prevents cells from entering mitosis when DNA is damaged, providing an opportunity for repair and inhibition of damaged cells proliferation. When mammalian cells contain damaged DNA, the p53 tumor suppressor and the Rb family of transcriptional repressors work together to downregulate a broad number of genes that encode proteins required for the G2 and M checkpoints. Elimination of these essential cell cycle proteins helps to keep the cells arrested in G2 [71,72]. Many studies have confirmed that miRNAs are involved in cell cycle arrest in the G2/M phase. In renal tubular cells, hypoxia induces G2/M arrest, leading to renal fibrosis via the miR-493-STMN-1 pathway [73]. miR-192 is responsible for the G2/M arrest in the proximal tubular epithelial cells after toxic injury by aristolochic acid [74]. On the other hand, evidence from endometriotic cells suggests that miR-210-3p attenuates the G2/M cell cycle checkpoint by inactivating the BRCA1 complex function in response to DNA damage under hypoxia via targeting the 3 untranslated region of BARD1 mRNA [75].
We also observed enrichment in the TBK1 gene set supporting the data of Barbie et al. who studied the epithelial lung cancer cell lines upon overexpression of an oncogenic form of the KRAS gene and knockdown of the TBK1 gene by RNAi. In experiments on human cancer cell lines, the lack of TBK1 induced apoptosis which was specifically dependent on oncogenic KRAS expression. In these cells, TBK1 activated NF-kappaB antiapoptotic signals involving c-Rel and BCL-XL, which were essential for survival and provided mechanistic insights into this synthetic lethal interaction. These observations indicate that TBK1 and NF-kappaB signaling is fundamental in KRAS-mutated tumors and establishes a general approach for the rational identification of codependent pathways in cancer [76].
The most important finding from this study was that high miR-154-5p and miR-154-3p expression levels were positively correlated with the EMT pathway in the HNSCC patients. It should be emphasized that the EMT pathway was enriched in all the localizations of HNSCC tumors. It had been confirmed before that miR-154 promoted EMT in prostate cancer cells. In several isogenic prostate cancers, the ARCaP and LNCaP lines that drive the EMT and bone metastasis in mice were significantly increased in populations with high miR-154 levels. Inhibition of miR-154-3p (miR-154*) in mice led to a decrease in bone metastases along with an increase in survival [77,78]. Interestingly, Chen et al. proved that overexpression of miR-154-5p in nasopharyngeal carcinoma (NPC) cell line models suppresses their invasion and migration. Their mechanistic studies showed the regulatory function of the said miRNA which, by directly binding the 3 UTR of kinesin-like protein 14 (KIF14) mRNA, inhibits NPC metastasis [62].
Based on the target prediction tools, we selected 24 genes that were potential targets for miR-154-5p, and the expression levels of those genes were negatively correlated with the expression levels of miR-154-5p in the HNSCC patients. In that group VPS4B, LRRC1, DMXL1, NPEPPS, XPO7, LIN54, CEACAM7, CNOT4, C12ORF29, WAC, LIN9, RNF138, GNAI3, and DCAF16 had not been described previously in the context of the HNSCC nor EMT process. Our results indicated that those genes were co-expressed, had physical interactions, and were involved in the ESCRT, cell cycle, as well as viral replication.
The connections with HNSCC and EMT have been reported for ten genes. Zhang et al. observed that MYO6 (myosin VI) was upregulated in OSCC and was associated with the regulation of cell cycle progression and apoptosis [79]. In gastric cancer, MYO6 acts as a pro-metastatic or pro-EMT gene and is directly regulated by miR-143, as well as by miR-145 [80]. Moreover, the identified potential target of miR-154-5p was OVOL1 (ovo-like transcriptional repressor 1) which was described as an important transcription factor in the mesenchymal-to-epithelial transition (MET) process [81,82]. It was proved that OVOL1 together with OVOL2 influenced the cellular phenotype by regulation of transcription factor ZEB1 and epithelial splicing regulatory protein 1 (ESRP1), which regulated mRNA splicing [82]. It has also been postulated that the family of transcription factors OVOL play a multiway function in cellular phenotype regulation by preventing EMT, driving MET, regulating the transition state (hybrid epithelial/mesenchymal), and maintaining a twostep process of transition EMT and MET [83]. We also observed that ABI1 (abl interactor 1) was negatively correlated with miR-154-5p. It has been indicated that low ABI1 expression was associated with metastasis and shorter survival of prostate cancer patients. Based on the in vitro studies, it has been revealed that ABI1 causes activation of the noncanonical WNT signaling and EMT pathways [84]. Another study showed that ABI1 has an inverse function, and its upregulation induced the EMT process and increased stem cell activity in breast cancer [85]. Furthermore, Fang et al. showed that the SOS1/EPS8/ABI1 complex was pivotal for ovarian cancer cells during the EMT process [86]. We also identified CDKN2B as a potential target. Reduced expression of CDKN2B was determined as a marker of disease recurrence in oral cancer [87]. In the case of hepatocellular carcinoma, reduced expression of CDKN1A and CDKN2B is connected with the EMT process, migration, and cell cycle progression [88]. Another identified target within this study was ESRP1. ESRP1 and ESRP2 play an important role in alternative splicing, and their expression is reduced during the EMT process [89]. It has been indicated that both were upregulated in dysplastic samples of OSCC, but their reduction was probably restricted in the cells that had a motile phenotype acquired during cancer invasion [90]. The next two miR-154-5p targets connected with EMT were CPEB3 and PLAGL2. Zhong et al. demonstrated that CPEB3 played a crucial role in the EMT process by regulation of the interaction between cancer cells and tumor-associated macrophages (TAMs). Its functions included inhibition of IL-6R/STAT3 signaling by binding to IL-6R mRNA in cancer cells as well as inhibition of the IL-6 molecule in TAMs by CPEB3. Finally, CPEB3 prevented secretion of CCL2 molecules in cancer cells and caused a reversion of the polarity of M2 macrophages [91]. Moreover, the role of CPEB3 in the progression of cervical cancer as well as in tumorigenesis and metastasis of lung adenocarcinoma was shown [92,93]. PLAGL2 (PLAG1-like zinc finger 2) seems to have an opposite role to CPEB3. It has been shown that PLAGL2 promoted EMT by USP37mediated deubiquitination of Snail1 in gastric cancer [94], as well as β-catenin-dependent regulation of ZEB1 in colorectal cancer [95]. Additionally, higher levels of PLAGL2 were associated with metastasis [95] and finally with worse patient survival [94]. It should be noted that Wang et al. identified PLAGL2 as well as MAPKAPK5-AS1 as the direct targets of miR-154-5p and the MAPKAPK5-AS1/PLAGL2/HIF-1α signaling loop that is responsible for progression and metastasis of hepatocellular carcinoma [96]. We indicated that the next potential target for miR-154-5p was PTBP3. Wu et al. observed that PTBP3 was responsible for migration through the regulation of E-cadherin, influencing the EMT process, and its higher expression was associated with shorter survival in NSCLC [97]. However, Hou et al. demonstrated that PTBP3 was associated with lymph node metastasis, advanced stages, and poor OS of breast cancer patients. PTBP3 directly regulates the EMT process and cancer stem cells by binding to the 3 UTR of ZEB1 [98]. The last two potential miR-154-5p targets were TMPRSS11D and NSD2. Both of them have been described in HNSCC. The first one, TMPRSS11D, was one of the nine hub genes significantly correlated with HNSCC, and its higher levels were associated with better survival [99]. The other one, NSD2, was identified as altered in 6% of HNSCCs and is one of the lysine methyltransferases [100]. It has been observed that the low expression level of NSD2 was correlated with reduced OS of HPV-positive patients [101]. Nevertheless, NSD2 was overexpressed in patients with advanced cancer, who had poorly differentiated tumors, and its levels were higher than in normal epithelia. NSD2 directly regulated the transcription of NIMA-related kinase-7 (NEK7), which is a cell cycle regulator [102]. Moreover, lower levels of NSD2 caused a higher expression of the E-cadherin protein and decreased N-cadherin as well as vimentin. NSD2 interacts with TWIST1 and causes upregulation of H3K36me2, which, in turn, is manifested with EMT promotion [103].
Finally, our analysis of immunological features in groups of patients with low or high expression of both miR-154 mature strands demonstrated some intriguing differences. The obtained results suggest that individuals with decreased miR-154-5p and miR-154-3p can be characterized by more favorable features. First of all, our assessment of the ESTIMATE score determined that these patients have lower tumor purity. Decreased values of this parameter were previously correlated with worse prognosis in, e.g., gastric cancers [104,105], colon cancer [106], and glioma [107]; however, it should be noted that the ESTIMATE score is calculated based on the stromal and immune scores. The large amount of stromal cells recruited from endogenous tissue promotes tumor growth and secrete many factors that influence angiogenesis, proliferation, and metastasis [108], though increased infiltration of immune cells is far more complicated to interpret.
Mandal et al. demonstrated that HNSCCs are among the most highly immuneinfiltrated cancer types with an especially abundant T cell fraction, which indicates a better prognosis and response to immunotherapeutic approaches [109]. It should be emphasized that the positive impact is associated with high fractions of CD3+, CD8+, and regulatory T cells [109,110]. Interestingly, overexpression of CD3 is strongly correlated with better clinical outcomes of specifically HPV-negative individuals [110,111]. The increased counts of CD8+ lymphocytes do not bear an unequivocal prognostic function when correlated with the HPV infection status [110,112]. Nevertheless, Balermpas et al. showed that high levels of both CD3 and CD8 cells were significantly associated with improved distant metastasisfree survival (DMFS), which could imply the presence of a systemic immunosurveillance mechanism suppressing the development of micrometastases [110]. It is worth mentioning that da Silva et al. discovered a link between CD8+ T cells and expression of the genes coding CTA [113]. Protein products of these sequences can be found in only one healthy tissue-in the testis, and in a variety of tumor types, where they induce an immunological response [114]. It has been proven that in HNSCC, high levels of CTAs associated with a good prognosis are correlated with elevated counts of CD8+ T cells [113]. Our results also indicate the link between better OS and an increased level of CD8+ lymphocytes in low miR-154 expression groups; however, the CTA scores are higher in the opposite groups of patients. We believe it is due to the fact that we worked on Thorssons' CTA score values, which were assessed using a pool of various CT genes, without distinguishing them into groups related to good or poor prognosis [36].
Analysis of immune cell fractions infiltrating the tumor showed a decreased number of macrophages in the groups of patients with low expression of both miR-154 strands. Kumar et al. proved that high density of TAMs, specifically of those with the M2-like phenotype, plays a crucial role in tumor progression, as well as metastasis, and correlates with worse clinicopathological features [115]. Interestingly, TAMs are characterized by downregulated expression of E-cadherin and elevated levels of mesenchymal markers, e.g., vimentin, Snail, and Slug, which could imply its involvement in the EMT [116]. A study by She et al. observed that CD163-positive macrophages (M2-like) play a pro-tumoral role and, when co-cultured with cancer cells, release TGF beta, epidermal growth factor (EGF) and upregulate ERK1/2, inducing tumor growth [116]. This type of macrophages inhibits proinflammatory M1-like ones through the production of, e.g., IL-10, TGF beta, as well as VEGF, leading to angiogenesis and tissue remodeling [117]. Additionally, IL-10 and TGF beta induce conversion of the M0 macrophages to the M2c type [118]. It should be emphasized that in HNSCC, pathways including TGF beta and FGFR are known to interact with EGFR signaling, which affects the wound healing process and could lead to EMT [119]. All of the above is reflected in our results, especially in the case of miR-154-5p. The high expression group had not only elevated macrophages, particularly M0 and M2, but also increased TGF beta response and elevated wound healing signatures, which confirms its substantial role in EMT and tumor progression. In addition, we determined that miR-154-5p reaches the highest expression level in the mesenchymal subtype of HNSCC in comparison to others, which seems to additionally corroborate the above statement. A recent study analyzing immune signatures in HNSCC proposed three different subtypes of these malignancies Immunity-H, -M, and -L (high, medium, and low, respectively). The first one was characterized by increased immune and stromal infiltration, low tumor purity, low stemness, as well as intratumor heterogeneity, genomic stability, and favorable prognosis in contrast to the immunity-L subtype [120]. These findings seem to support our obtained results. Yoshihara et al. suggested that in HNSCC, tumor purity is correlated with mutation rates, which is especially interesting while taking into account higher silent and non-silent mutation rates discovered by us in a group of individuals with decreased levels of miR-154-5p [47]. They also implied that in some low-purity neoplasms, reduction in the T > A substitutions indicates an impact of the tumor microenvironment on mutational processes or that mutation types can alter the stromal and immune infiltrations in tumors [47]. Furthermore, the above study suggested that the stromal expression pattern could overlap with the mesenchymal phenotype of tumors, which, in addition to the strong correlation with tumor purity, could result in incorrect classification of an increased presence of tumorassociated stroma as an EMT process [47]. Nonetheless, our knowledge of the miR-154-5p involvement in the abovementioned transition, along with the data obtained during this research, excludes the possibility of the aforementioned misinterpretation.
The Immunity-H subtype is associated with a higher rate of HPV infections in comparison to the Immunity-L subtype (40% vs. 22%) [120]. It has been proven that HPV-positive patients with oropharyngeal squamous cell carcinoma (OPSCC) who are current smokers or have tobacco use history tend to have a dysregulated ncRNA expression landscape, including downregulation of miR-154-5p, and unfavorable treatment outcomes [121]. A study by Huang et al. indicated that smokers and lung cancer patients also manifest with low levels of miR-154 mature strands. Significant differences between these two groups of individuals and controls suggest the miR-154-5p diagnostic potential and imply an important role in the biology of cigarette smoke-induced lung cancer [122]. On the other hand, analysis of the miR-154-3p levels in plasma samples in groups of healthy smokers, individuals with lung granuloma and lung adenocarcinoma showed its overexpression in both types of patients. Those significant differences suggest that this molecule could help discriminate adenocarcinoma from lung granuloma in the future [123]. Interestingly, HNSCC patients with high molecular smoking signatures have an elevated number of mutations and decreased T cell infiltration, which is associated with local immunosuppression, low levels of cytotoxicity, and poor survival [109].
The broad network of miR-154 interactions includes many molecular pathways essential for tumor development and progression. Further research to explore the role of this molecule, as well as its diagnostic and therapeutic potential, could yield results relevant to the quality and life expectancy of HNSCC patients.
In conclusion, our results indicate that miR-154 plays an essential role in HNSCC development and progression, presumably through its influence on the Epithelial-to-Mesenchymal transition process and the regulation of phenotype of both cancer and immune cells. miR-154-5p targets the genes involved in the ESCRT, cell cycle, and viral replication, which are connected with the EMT process in other cancers as well. miR-154 has a potential to be used as a biomarker describing clinical, molecular, and immunological features of HNSCC patients. Moreover, the implementation of miR-154 in diagnosis and therapy might lead to the progress in the discovery of modern neoplasm treatment methods. However, it should be noted that the above data were obtained from in silico analyses and require further validation. miR-154 is a molecule with a tremendous medical potential whose exploration may bring significant improvement in understanding HNSCC, and thus more effective treatment development.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biomedicines9121894/s1, Table S1: List of processes based on the GSEA analysis of MSigDB gene sets enriched in HNSCC patients with lower and higher levels of miR-154-5p in all localization as well as depending on the localization in oral cavity, larynx or pharynx, Table S2: Correlation between miR-154-5p and 24 genes selected based on miRNA targets with score, function and interactions between targets using the GeneMANIA tool.

Institutional Review Board Statement:
The study was based on the analysis of freely available datasets and did not need any ethics committee's approval and does not violate the rights of other persons or institutions.
Informed Consent Statement: Not applicable. This study used the data from TCGA, which are freely available and do not violate any rights of people or institutions.
Data Availability Statement: All the data are available online with common access. The data analyzed during this study are available from the corresponding author on reasonable request.