Comprehensive Landscape of Prognostic Significance and Immune Characteristics of Myosins in Squamous Cell Carcinoma of the Head and Neck

Myosin superfamily, a large and diverse family of molecular motors important for cell motility and migration, has been illustrated to play contradictory roles during the development of several kinds of tumors. However, the function and prognostic values of MYOs in head and neck squamous cell carcinoma (HNSCC) still remain largely unknown. In the current manuscript, the expression levels and clinical data of MYOs in HNSCC were investigated by online databases, including Oncomine, GEPIA, GEO, TCGA, HPA, UALCAN, Kaplan-Meier plotter, and CancerSEA; we found that the expression levels of MYO1B, MYO5A, and MYO10 were significantly elevated in HNSCC tissues, which were also correlated with the unfavorable overall survival (OS) of the patients. Furthermore, MYO1B/MYO5A/MYO10 interacting genes were identified, and the protein-protein interaction (PPI) networks were constructed by STRING and GeneMANIA. The enrichment analysis revealed that MYO1B/MYO5A/MYO10 associated genes mainly participated in cell metastasis and EMT processes, which were also confirmed by cell functional experiments. At last, the ssGSEA method was conducted to investigate the extent of immune cell infiltration, and we found that both the expression of MYO1B/MYO5A/MYO10 were closely correlated with the infiltration of immune cells in HNSCC. These findings implied that MYO1B, MYO5A, and MYO10 as novel prognostic factors for HNSCC and provided new strategy for HNSCC treatment.


Introduction
As the most frequent type of head and neck cancer and the seventh most common cancer worldwide, head and neck squamous cell carcinoma (HNSCC) presents extremely high morbidity and mortality and contributes to severe healthy burden [1]. HNSCC includes malignancies in the regions of the oral cavity, oropharynx, nasopharynx, hypopharynx, and larynx [2], and the occurrence of HNSCC is often thought to be correlated with HPV virus infection [3]. Regional neck lymph node metastasis is a prone feature of HNSCC; once evolved into distant metastasis, it will seriously affect patient survival and prognosis [4]. Therefore, there is an urgent need to develop prognostic biomarkers and to identify new therapeutic targets for HNSCC.
Recent studies revealed that various MYOs played crucial roles in tumorigenesis and cancer development [8]. For instance, MYO1B was illustrated to contribute to cell proliferation, migration, and invasion and enhanced the activities of MMP1/MMP9 in cervical cancer [9]. Knockdown the expression of MYO10 alleviated tumor invasion and metabolic stress responses in glioblastoma [10]. MYO5B expression was reported to be upregulated in pheochromocytoma and paraganglioma tissues, and MYO5B mutation (p.L587P, p.G1611S, and p.R1641C) was demonstrated to be responsible for cancer cell proliferation and migration [11]. In HNSCC, increased expression of MYO1B was shown to aggravate cell migration and lymph node metastasis by enhancing cell motility [12], and silencing of MYO6 contributed to the inhibition of cell proliferation via regulation of cell cycle and apoptosis in oral squamous cells [13]. However, the expression and prognostic values of other MYOs in HNSCC have not been comprehensively demonstrated.
To the best of our knowledge, bioinformatics analysis has yet to be applied to explore the roles of MYOs in HNSCC. The present study assessed the correlation between MYO expression and the prognostic indicators of HNSCC by using various online analysis databases and identified the potential biomarkers for the treatment of HNSCC.

UALCAN.
To investigate the expression of MYOs in cancer stages and nodal metastasis status of HNSCC, UAL-CAN database (http://ualcan.path.uab.edu/) [19] which was designed to analyze the expression level of genes in TCGA and the clinical data of patients was applied.

Kaplan-Meier Plotter Database and Receiver Operating
Characteristic (ROC) Analysis. The prognostic significances of MYO expression in terms of overall survival (OS) information in HNSCC were demonstrated by utilizing Kaplan-Meier plotter database (http://kmplot.com/analysis/) through dividing patients to high or low expression groups of MYOs by auto selected best cut-off option. Diagnostic values of MYOs were elucidated by ROC analysis by using the pROC package based on clinicopathological parameters from TCGA, and area under the curve ðAUCÞ > 0:8 was considered as an ideal biomarker for distinguishment.
2.6. Correlation Analysis and CancerSEA Database. Correlation between the mRNA expression of MYOs was detected by using Pearson's correlation coefficient method and Corrplot package. CancerSEA (http://biocc.hrbmu.edu.cn/ CancerSEA/) was applied for illustration of the functional states of cancer cells at the single-cell resolution [22].

Functional Enrichment Analysis.
For functional enrichment analysis, we identified differentially expressed genes (DEGs) in MYO-high and MYO-low groups (|log 2FC | >1 and adjust p value < 0.05 as thresholds) based on TCGA data, and the top 150 positively associated DEGs were identified, then subjected to Gene Ontology (GO) enrichment included biological process (BP), molecular function (MF), and cellular component (CC) enrichment, as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment with clusterProfiler package. For gene set enrichment analysis (GSEA) analysis, Hallmark gene set (h.all.v7.2.symbols.gmt) from MSigDB database (https://www.gsea-msigdb.org) was utilized with the clusterProfiler package visualized by ggplot2 package. Adjust p value < 0.05, FDR < 0:25, and |NES | >1 were considered significant enrichment.
2.9. Single-Sample Gene Set Enrichment Analysis (ssGSEA). The immune cell infiltration levels were measured by the ssGSEA method using GSVA package (http://www.bioconductor.org/ packages/release/bioc/html/GSVA.html) as described [24]. Significance was determined by the Wilcoxon rank sum test and Spearman correlation method.

Statistical Analysis.
For bioinformatic analysis, statistical analyses were carried out using the R package. Student's t-test and one-way ANOVA together with Tukey Kramer post hoc testing were used in the cellular functional experiments.

Transcription Levels of MYOs in HNSCC.
In order to investigate the roles of MYOs in HNSCC, Oncomine data-base was selected to analyze 25 MYO expressions in cancers. As shown in Figure 1 (Figure 2(a)), as well as the expression data in GSE31056 dataset of GEO database (Figure 2(b)). All these results confirmed the aberrant transcription levels of MYO1B, MYO5A, MYO5C, and MYO10. Furthermore, the protein levels of MYO1B, MYO5A, MYO5C, and MYO10 in HNSCC tumor tissues and normal tissues were examined by Human Protein Atlas (HPA) database; consistently, the protein levels of MYO1B and MYO10 were dramatically enhanced in tumor tissues compared with normal tissues (Figure 2(c)).

Correlation and Interaction Network of MYOs in HNSCC.
Moreover, we examined the correlation between MYOs by the Pearson correlation method, and we found that MYO1B showed a positively correlation with MYO5A and MYO10, but negatively correlated with MYO5C (Figure 4(a)). In addition, the correlation between MYO5C and MYO5A or MYO10 was not significant. Taken together, the aforementioned details of MYO5C on expression, survival data, and diagnostic value, we considered that it may not be sufficient as a stable indicator of HNSCC; therefore, we mainly illustrate the effects and values of MYO1B, MYO5A, and MYO10 in the following studies. Due to the correlation between MYO1B, MYO5A, and MYO10, we constructed the protein-protein interaction network of the three MYOs, and 50 associated proteins in STRING database (Figure 4(b)) and 20 associated proteins in GeneMANIA database were identified (Figure 4(c)).

Functional States of MYOs in CancerSEA Database.
After certification of MYO expression, we further explored the potential mechanisms and related biological functional processes of MYO1B, MYO5A, and MYO10. As shown in Figure 5(a), MYO expression distributions in single-cell resolution were analyzed by CancerSEA, and the cells with high expression of MYO1B/MYO5A/MYO10 were tended to cluster together, which implied the contribution of MYOs in malignant progression. Furthermore, we found that MYO1B, MYO5A, and MYO10 were positively correlated with functional states including metastasis, invasion, and EMT ( Figures 5(b) and 5(c)).

Functional Enrichment Analysis of MYOs Associated
Genes. In order to evaluate the functional processes and pathways of MYOs, we identified differentially expressed genes (DEGs) in high and low MYO expression groups (criteria set of |log 2FC | >1 and adjust p value < 0.05), and the top 150 positively associated DEGs were selected to perform GO (Figure 6(a)), KEGG (Figure 6(b)), and GSEA enrichment analyses (Figure 7). The involvement of MYO1B associated DEGs was mainly in the receptor and integrin binding, focal adhesion and cell-substrate junction, epithelial cell differentiation, matrix organization, etc. The enrichment of MYO5A associated DEGs was in receptor binding and peptide activity, cell junction and collagen-containing matrix, peptide cross-linking, and epidermal cell differentiation. GO enrichment of MYO10 associated DEGs was mainly involved in molecular binding and matrix constituent, focal adhesion and collagen-containing matrix, cell adhesion, and matrix organization (Figure 6(a)). In KEGG analysis results, we noticed that MYO1B, MYO5A, and MYO10 associated DEGs were enriched in PI3K-Akt signaling, focal adhesion process, and ECM-receptor interaction which were all related to tumor metastasis [25][26][27]. Moreover, GSEA enrichment data suggested that MYO1B (Figure 7(a)),  (Figure 7(b)), and MYO10 (Figure 7(c)) were enriched in the pathways of EMT, apical junction, and myogenesis.

Silencing of MYOs Inhibits Cell Migration, Invasion, and EMT.
To verify the effect of MYOs in metastasis and EMT process, we used siRNAs specific targeted MYOs to knockdown MYO expression, and the efficiencies of siRNAs were confirmed (Figure 8(a)). Transwell assay results showed that silencing of MYO1B, MYO5A, and MYO10 significantly inhibited FaDu cell migration and invasion (Figures 8(b) and 8(c)). Furthermore, we examined the expression of EMT markers including E-cadherin, N-cadherin, and Vimentin in MYOs silenced FaDu cells, and we found that siRNA treatment greatly increased E-cadherin expression and downregulated the protein levels of N-cadherin and Vimentin, which meant that silencing of MYO1B, MYO5A, and MYO10 both alleviated the EMT potential of FaDu cells (Figure 8(d)).

Association between MYOs and Immune Cell Infiltration in HNSCC.
Immune infiltration and tumor immune microenvironment have been demonstrated to play essential roles in HNSCC [28,29]. To estimate the association between MYO expression and immune cell infiltration, single sample gene set enrichment analysis (ssGSEA) was conducted. As shown in Figure 9(a), MYO1B, MYO5A, and MYO10 were positively associated with neutrophils, T gamma delta cells (Tgd), and T central memory cell (Tcm) cell infiltration, whereas negatively correlated with pDC, CD8 T cells, and cytotoxic cells. The infiltration of neutrophils was significantly aggravated in the MYO1B/MYO5A/MYO10 high expression group compared to the MYO-low expression group in TCGA database (Figure 9(b)). The correlation between MYOs and neutrophils was also shown in scatter plots (Figure 9(c)).

Discussion
Despite advances in diagnostic detection and surgical techniques, the high rates of recurrence and metastasis are still considered as mainly restricting factors in HNSCC treatment. Therefore, in-depth exploration of the crucial mechanisms and target molecules related to the development and metastasis of HNSCC is of great significance for the treatment of HNSCC. In the current study, for the first time, we comprehensively analyzed the expression details and prognostic significances of MYOs and suggested MYO1B, MYO5A, and MYO10 as effective clinical biomarkers and potential medical targets for HNSCC.
Although members of the myosin superfamily are involved in almost all aspects of human life [30], the roles of MYOs in cancers especially in HNSCC still remain to be elucidated. Benefiting from rapid advances in bioinformatics and sequencing technologies, online information databases have been significantly expanded, providing us with great help to analyze and study potential molecular markers. Previous researches revealed that MYO1B expression levels were elevated in cervical cancer [9], prostate cancer [31], colorectal cancer [32], and oral tongue cancer [12]. MYO5A was also found to be increased in glioblastoma [33] and esophageal carcinoma [34] tissues compared to the normal tissues. Aberrant levels of MYO10 were observed in breast cancer [35] as well as squamous cell carcinoma of the lung [36]. In the current study, with the utilization of databases including Oncomine, GEPIA, and GEO (GSE31056 dataset), we comprehensively analyzed 25 MYO transcription levels and obtained 4 intersection genes (MYO1B, MYO5A, MYO5C, MYO10) among these three databases. Consistent with previous reports about the elevated expression in other cancers, MYO1B, MYO5A, and MYO10 levels were increased in HNSCC tissues compared with adjacent normal tissues whereas MYO5C expression was found to be downregulated. Moreover, high expression of MYO1B, MYO5A, and MYO10 was shown to associate with the unfavorable overall survival of HNSCC patients, while the prognostic value of MYO5C was not significant as well as the moderate ROC diagnostic value. These results suggested the great potential of MYO1B, MYO5A, and MYO10 as prognostic biomarkers for HNSCC.
Recently, MYO1B was reported to aggravate colorectal cancer metastasis via enhancing rearrangement of F-actin and focal adhesion assembly mainly through targeting RhoA [32]. MYO5A was shown to be regulated by snail and contributed to cancer cell migration and invasion [37]. Cao et al. revealed that MYO10 aggravated the aggressiveness and metastasis of breast cancer cells through invadopodial formation [38]. To further explore the underlying mechanisms of MYOs in HNSCC, we used CancerSEA database to investigate MYO expression and correlated functional states. We found that the biological processes of metastasis, invasion, and EMT were extremely correlated with the levels of MYO1B, MYO5A, and MYO10. We next divided TCGA patients into high and low MYO expression groups and identified the most associated DEGs, followed by functional enrichment analysis. By synthesizing the findings of GO, KEGG, and GSEA analyses, we found similarities in the enrichment results of MYOs associated DEGs. For instance, cancer development-related PI3K-Akt signaling, focal adhesion signaling, integrin-ECM pathways, and EMT processes were all enriched in MYO1B/MYO5A/MYO10 associated DEGs. To strengthen the conclusion, we also performed cellular functional experiments, and we found that silencing the expression of MYO1B, MYO5A, and MYO10 alleviated the abilities of cell migration and invasion, and protein levels of EMT markers including E-cadherin, N-cadherin, and Vimentin were all regulated by MYO-siRNA administration.
Although the expression and functions of several MYOs in tumors have been discussed, the effects of MYOs on immune cell infiltration still remain largely unknown. Accumulating evidences suggested that immune infiltration and tumor microenvironment were extremely important in tumorigenesis, tumor development, and metastasis [39,40]. In the current manuscript, we evaluated the association of MYO1B, MYO5A, and MYO10 expressions and immune cell infiltration, and we found that various immune cells such as neutrophils, T gamma delta cells (Tgd), T central memory cells (Tcm), pDC, and CD8 T cells were positively/negatively correlated with MYO1B, MYO5A, and    14 Journal of Immunology Research MYO10 expressions, which meant that MYOs may participate in the immune processes. It is worth noting that among the immune cells, CD8 T cells and cytotoxic cells were significant negatively correlated with expression of MYO1B, MYO5A, and MYO10, suggesting that high expression of MYOs may inhibit the effect of antitumor immunity. Most importantly, we noticed that neutrophil infiltration exhibited great correlation with both expression of MYO1B/ MYO5A/MYO10. Recently, neutrophils were found to release chromatin DNA filaments coated with granule proteins to form neutrophil extracellular traps (NETs) and aggravated tumor development and metastasis [41][42][43], which suggested that MYO1B, MYO5A, and MYO10 may associated with NET formation therefore affected EMT process and tumor metastasis.

Conclusions
In conclusion, this manuscript provided comprehensive analyses about the expression and function of MYOs in HNSCC, suggested MYO1B, MYO5A, and MYO10 as potential prognostic biomarkers in HNSCC, and revealed the potential novel roles of MYO1B, MYO5A, and MYO10 in tumor immune microenvironment.

Data Availability
Publicly available datasets were analyzed in this study. The data used and/or analyzed during the current study are available from the corresponding author on reasonable request.