MicroRNA-Based Markers of Oral Tongue Squamous Cell Carcinoma and Buccal Squamous Cell Carcinoma: A Systems Biology Approach

Objective Oral tongue squamous cell carcinoma (OTSCC) and buccal squamous cell carcinoma (BSCC) are the first and second leading causes of oral cancer, respectively. OTSCC and BSCC are associated with poor prognosis in patients with oral cancer. Thus, we aimed to indicate signaling pathways, Gene Ontology terms, and prognostic markers mediating the malignant transformation of the normal oral tissue to OTSCC and BSCC. Methods The dataset GSE168227 was downloaded and reanalyzed from the GEO database. Orthogonal partial least square (OPLS) analysis identified common differentially expressed miRNAs (DEMs) in OTSCC and BSCC compared to their adjacent normal mucosa. Next, validated targets of DEMs were identified using the TarBase web server. With the use of the STRING database, a protein interaction map (PIM) was created. Using the Cytoscape program, hub genes and clusters within the PIM were shown. Next, gene-set enrichment analysis was carried out using the g:Profiler tool. Using the GEPIA2 web tool, analyses of gene expression and survival analysis were also performed. Results Two DEMs, including has-miR-136 and has-miR-377, were common in OTSCC and BSCC (p value <0.01; |Log2 FC| > 1). A total of 976 targets were indicated for common DEMs. PIM included 96 hubs, and the upregulation of EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, and HSPA5 was significantly associated with a poor prognosis in the head and neck squamous cell carcinoma (HNSCC), while NTRK2, HNRNPH1, DDX17, and WDR82 overexpression was significantly linked to favorable prognosis in the patients with HNSCC. “Clathrin-mediated endocytosis” was considerably dysregulated in OTSCC and BSCC. Conclusion The present study suggests that has-miR-136 and has-miR-377 are underexpressed in OTSCC and BSCC than in normal oral mucosa. Moreover, EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, HSPA5, NTRK2, HNRNPH1, DDX17, and WDR82 demonstrated prognostic markers in HNSCC. These findings may benefit the prognosis and management of individuals with OTSCC/BSCC. However, additional experimental verification is required.


Introduction
Oral cancers (OCs) are a frequent type of malignancy in the head and neck region with a poor treatment outcome. Based on a previous report, the GLOBOCAN study estimated that 377,713 of the world's population would be afected by lip and oral cavity cancer in 2020, leading to 177,757 deaths [1]. Te main risk factors for OCs are the use of tobacco and alcohol, exposure to UV radiation, and infection with the human papillomavirus (HPV) and Epstein-Barr virus (EBV) [2]. In addition, the aggressive nature of OC cells is strongly correlated with matrix metalloproteinase (MMP)-2, -9, and -13 [3]. Radical resection, radiotherapy, and chemotherapy are treatment approaches for OC patients. Among all the OC cases, 90% are classifed as squamous cell carcinoma (SCC) [4]. Oral tongue SCC (OTSCC) is the most frequent subtype of OC and has the worst prognosis among all oral squamous cell carcinomas (OSCCs) [5][6][7][8]. Besides, buccal squamous cell carcinoma (BSCC) was reported as the second most common tumor in the oral cavity [9]. Despite advancements in therapeutic methods, the 5-year survival rate of OSCC patients has not improved efciently. When the illness is discovered at stage IV, it drops from 80% to 30% when compared to primary OSCC [10]. By the appearance of clinical signs, up to 50% of OSCC patients are identifed after they are already advanced [11]. In order to battle OSCC and increase patient survival rates, it is necessary to suggest more efective treatment options [12].
MicroRNAs (miRNAs) are noncoding and small RNAs (20-23 nucleotides) that regulate the transcription of their target genes; these upstream regulators bind to their complementary sequences at 3′ UTR of their mRNA targets, leading to mRNA degradation or translation inhibition [13][14][15]. New research revealed that miRNAs might drastically dysregulate the expression of genes linked to ER stress, necroptosis, pro-and antiapoptotic activity, and oncogenes [16,17]. Terefore, miRNAs showed regulatory roles in critical biological procedures, including the cell cycle process, apoptosis and necroptosis, proliferation, and diferentiation [18]. Consequently, miRNAs were defned as valuable biomarkers for the prognosis and diagnosis of patients with cancer, which have demonstrated satisfying results when assigned as therapeutic targets in cancer [12].
Accumulating evidence suggests that miRNAs play a critical role in the pathogenesis and prognosis of patients with head and neck squamous cell carcinoma (HNSCC). An earlier study by Dioguardi et al. [19] found a strong correlation between miR-31 overexpression and a poor prognosis in HNSCC patients. Additionally, the miR-196 family [20] and miR-155 [21] have been identifed as potential predictive indicators of survival for HNSCC. Besides, Dioguardi et al. [22] showed the underexpression of miR-195 in patients with HNSCC and its excellent potential as an independent prognostic survival marker for HNSCC.
Te present study hypothesized that miRNAs might participate in the tumorigenesis of OTSCC and BSCC. It was suggested that there are common diferentially expressed miRNAs (DEMs) in the OTSCC and BSCC tissues compared to the normal oral mucosa, leading to abnormal expression of their target genes. Tese targets may be involved in several signaling pathways and biological processes (BPs), which may help elucidate the mechanisms underlying OSCC. Some other target genes might act as prognostic markers in patients with OSCC. Tus, utilizing the orthogonal partial least squares (OPLSs), common DEMs in OTSCC and BSCC in comparison to their respective healthy tissues were revealed. Ten, validated targets of common DEMs were identifed, a protein-protein interaction (PPI) network was constructed, hub genes within the network were illustrated, and the prognostic roles of hubs in HNSCC were evaluated using Kaplan-Meier curves. Furthermore, the enriched pathways and BP terms associated with the main clusters in the PPI network were identifed.
Te dataset GSE168227 from the Gene Expression Omnibus (GEO), available at https://www.ncbi.nlm.nih.gov/ geo/, was considered for reanalyzing to examine the present hypothesis. Te dataset was developed by Rajan et al. [23] to monitor the miRNA expression pattern in OTSCC, BSCC, and their adjacent oral mucosa to achieve a novel prognostic miRNA signature for OSCC. All tissue samples were from patients at the Regional Cancer Center's Head & Neck Clinic (Tiruvananthapuram, India). Patients with any chronic systemic condition or those who had already had cancer therapy were not allowed to participate in the trial. Control samples were attained from the patients undergoing oral and maxillofacial surgery for noncancer illnesses. All tissue specimens were immediately snap-frozen and stored in liquid nitrogen.

Dataset Recovery and Statistical Analysis.
Te miRNA expression dataset GSE168227 [23], based on the platform GPL8227 (Agilent-019118 Human miRNA Microarray 2.0 G4470B), was downloaded as a TXT fle from Gene Expression Omnibus (GEO), available at https://www.ncbi. nlm.nih.gov/geo [24]. Te dataset contained 48 oral tissue samples, including OTSCC (n � 16), BSCC (n � 14), normal tongue samples (n � 8), and normal buccal specimens (n � 10). For feature selection, R programming (version 4.0.2) was used. DEMs in OTSCC and BSCC were found using OPLSs in comparison to the comparable normal tissues. A cut-of condition was set to p value <0.01 and | Log2 fold change (FC)| > 1 [25,26]. Common DEMs involved in the malignant transformation of normal oral mucosa to OTSCC and BSCC were indicated. Furthermore, the Shiny apps web-based tool [27] showed the volcano plot of miRNAs in the dataset GSE168227.

PPI Network Analysis Based on Common DEM Targets.
Te experimentally validated targets of common DEMs were identifed using the TarBase version 8 database, available at https://dianalab.e-ce.uth.gr/html/diana/web/index.php?r= tarbasev8/index [28]. Possible interactions among the targets with a combined score of 0.4 [29] were indicated using the STRING knowledge tool, available at http://string-db. org/ [30]. After removing disconnected nodes inside the protein interaction map (PIM) [31], the connected network was downloaded as a TSV format and subsequently imported into the Cytoscape 3.9.1 software [32], available at https://www.cytoscape.org to perform structural analysis. Hub nodes were awarded to proteins whose degree and betweenness were above two times the average and whose closeness exceeded the mean of the PPI network's nodes. In line with our earlier research [33], utilizing the MCODE plugin, the primary clusters within the PIM were also highlighted.

Gene Set Enrichment
Analysis. Signifcant molecular pathways and Gene Ontology (GO) terms involved in the malignant transformation of normal oral mucosa to OTSCC and BSCC were unraveled using g:Profler tool, available at https://biit.cs.ut.ee/gprofler/gost [34]. Te cut-of condition was set to the false discover rate (FDR) < 0.05 and the number of entities ≥10, following our previous study [35]. Te results from two primary pathway sources, including Reactome [36] and the Kyoto Encyclopedia of Genes and Genomes (KEGGs) [37] databases, were considered for pathway enrichment analysis. Notably, common DEMs targets were taken into account for enrichment analyses of cellular component (CC) and molecular function (MF). At the same time, pathway and BP annotation enrichment analyses were given to the genes inside the clusters as input data [15].

Prognostic Role of Hub Genes.
Regarding the critical role of hub genes in the pathogenesis of OTSCC and BSCC, the prognostic impact of hubs in the HNSCC was evaluated using the Gene Expression Profling Interactive Analysis 2 (GEPIA2) database [38], available at http://gepia2.cancerpku.cn/#index. By reanalyzing the raw RNA-seq expression data of malignant and healthy tissue samples from the TCGA and GTEx datasets, the GEPIA2 generates Kaplan-Meier curves. Prognostic indicators were defned as genes with the log-rank test and the hazard ratio (HR) p values < 0.05. Te prognostic role of the combination of genes was also evaluated [39].

Gene Expression Evaluation of Prognostic Markers.
Te mRNA expression patterns of prognostic genes were evaluated in HNSCC tissues (n � 519) and healthy samples (n � 44). It was performed by boxplot analysis provided by the GEPIA2 database [38].

DEMs in OTSCC and BSCC.
Two prediction models were created using OPLSs to fnd DEMs between the OTSCC and BSCC tissues in comparison to their healthy counterparts. Each model was built using 369 variables and 24 samples. At p value <0.01 and Log2 FC > 1, a total of 10 DEMs were indicated in OTSCC compared to the healthy oral mucosa (R2X � 0.494; R2Y � 0.71; and Q2 � 0.264). Besides, seven DEMs were identifed in BSCC than in control samples (R2X � 0.526; R2Y � 0.715; and Q2 � 0.411) (Figure 1(a)). All DEMs in OTSCC and BSCC are listed in Table 1. Two DEMs, including has-miR-136 and has-miR-377, were common in two subtypes of OSCC. Volcano plots showed DEMs in the studies groups based on −Log 10 p value and Log 2 FC (Figure 1(b)).

Hub Genes, Modules, Pathways, and GO terms. Te
TarBase database detected 976 genes as experimentally validated targets for common DEMs. Te list of targets was used as input data in the STRING database to construct a PPI network. Single nodes were removed from the PIM, and subsequently, a connected network including 932 proteins and 6682 interactions was imported into Cytoscape. Te average degree, betweenness, and closeness value was calculated as 14.33, 0.0023, and 0.33, respectively. Ninety-six nodes were then identifed as hub genes linked to the pathogenesis of OTSCC and BSCC. (  (Figure 2). Te most important pathways implicated in the carcinogenesis of OTSCC and BSCC also included "pathway in cancer" (KEGG:05200), "proteoglycan in cancer" (KEGG:05205), "bladder cancer" (KEGG:05219), and "clathrin-mediated endocytosis" (REAC: R-HSA-8856828). Te two most essential BPs promoting malignant transformation in the tongue and buccal area were "cell death" (GO:0008219) and "apoptotic process" (GO: 0006915). Moreover, "neoplasm" (GO:0005654) and "nuclear lumen" (GO:0031981) CCs were considerably afected in OTSCC and BSCC. Moreover, "Protein-containing complex binding" (GO:0044877) and "RNA binding" (GO:0003723) were the most enriched terms in the category of MFs. Figure 3 demonstrates the most signifcant pathways and GO terms dysregulated in the pathogenesis of OTSCC and BSCC.

Prognostic Markers.
Kaplan-Meier curves showed that the overexpression of EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, and HSPA5 was signifcantly associated with a dismal prognosis in the patients with HNSCC. EIF2S1 showed the most negative marker with the criteria of HR � 1.7, log-rank test p value � 0.00016, and HR p value � 0.00019. Additionally, an improved prognosis in HNSCC was associated with the overexpression of NTRK2, HNRNPH1, DDX17, and WDR82. With an HR, log-rank p value, and HR p value of 0.71, 0.011, and 0.011, respectively, NTRK2 was the most positive marker. Terefore, these markers might be assigned as drug targets in patients with OTSCC and BSCC. Figure 4 presents the survival analysis of prognostic markers. Figure 5 shows interactions between hub genes. Te combination of EIF2S1, CAV1, RAN, and ANXA5 showed a considerable negative signature with the criteria of HR � 1.6, log-rank p value � 0.00028, and HR p value � 0.00032 (Table 3).

Discussion
Te incidence of OTSCC and BSCC is considered the frst and second highest among all oral cancers, respectively [40,41]. MMP-13 inhibitors [42] may be efective treatments for improving survival rates in patients with OTSCC and BSCC because of the substantial role MMP-13 plays in the invasion of OC cells. However, unraveling the most critical genes, molecular pathways, and GO annotations mediating the malignant transformation of normal oral mucosa to OTSCC and BSCC might be helpful in treating OSCC. Gene set enrichment analysis showed that the "Clathrin-mediated endocytosis" pathway (REAC: R-HSA-8856828) signifcantly mediates the malignant transformation of normal oral mucosa to OTSCC and BSCC. "Clathrin-mediated endocytosis" (CME) is an endocytic process that regulates the expression of plasma membrane receptors and infuences the pathways that lead to those receptors' downstream signals [43][44][45]. Xiao et al. [46] showed that ERK1/2 phosphorylates the FCH/F-BAR and SH3 domaincontaining protein (FCHSD2), leading to enhanced clathrin-coated pit (CCP) initiation and CME, resulting in decreased cell-surface EGFR expression and reduced proliferation and migration of nonsmall cell lung cancer cells. Te present study illustrated that hsa-miR-377 and hsa-miR-136 are commonly downregulated in OTSCC and BSCC compared to the normal oral mucosa. Sun et al. [47] reported a signifcant downregulation of miR-377-3p in nonsmall cell lung cancer (NSCLC) tissues compared to their corresponding healthy lung specimens. Besides, Sun et al. [47] demonstrated a remarkable negative correlation between the expression of oncogene E2F3 mRNA and miR-377-3p in lung tissues using Pearson correlation analysis (r2 � 0.3614, p < 0.0001), suggesting the tumor suppressive role of miR-377-3P by downregulating the E2F3. Te authors noted that the elevation of the long noncoding RNA (lncRNA) NEAT1 prevented the increased cell death caused by miR-377-3p. NEAT1 overexpression was linked to carcinogenesis and cancer development, and NEAT1 has been described as an oncogenic gene in several malignancies [48,49]. Evidence suggests that estrogen receptor-α36 (ERα36) plays a signifcant role in the tumorigenesis of luminal subtypes of breast cancer. In this regard, the enhanced expression of ERα36 is associated with disease development and drug resistance [50][51][52]. Tiebaut et al. [53] showed a negative correlation between the expression of the miR-136-5p and ERα36 in breast cancer cells. Te miR136-5p mimic transfection in MCF-7 cells diminished the ERα36 expression.
OSCC is the most common HNSCC [54], counting for 95% of all HNSCC cases [55]. Furthermore, the GEPIA2 (or GEPIA) database for HNSCC is commonly used to explore the prognostic efect of genes in OSCC and/or to evaluate the mRNA expression levels of genes in OSCC  compared to healthy controls. Using immunohistochemical analysis, Sun et al. [56] compared the protein expression levels of AKT1 and PLK1 in OSCC tissues to normal oral specimens. Subsequently, Sun et al. [56] used the GEPIA web server to validate their experimental results; this included evaluating the mRNA expression levels of AKT1 and PLK1 in HNSCC than in healthy tissues and analyzing the prognostic efect of the genes in patients with HNSCC. Fang et al. [57] identifed several genes with a correlation score >0.8 in a PPI network associated with the pathogenesis of OSCC. Next, the authors used the GEPIA tool to evaluate the prognostic role of the genes in HNSCC. Furthermore, Dai et al. [58] performed a weighted gene comethylation network analysis (WGCNA) to identify hub modules and CpG sites correlated with OSCC. Ten, Dai et al. [58] used the GEPIA for HNSCC to conduct a Kaplan-Meier survival analysis to investigate the possible predictive signifcance of the numerous hub CpG site-associated genes. Te GEPIA2 for HNSCC was used in the present study to investigate the potential prognostic impact of the hub genes in OTSCC and BSCC. Gene expressions of the prognostic markers were also evaluated using GEPIA2 for HNSCC. Based on the targets of hsa-miR-377 and hsa-miR-136, a PPI network of 932 genes and 6682 edges was created here. EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, and HSPA5 overexpression was substantially related to a poor prognosis in patients with HNSCC, and 96 hub genes showed a striking centrality within the PPI network. Besides, the overexpression of NTRK2, HNRNPH1, DDX17, and WDR82 was linked to a favorable prognosis in HNSCC patients (the log-rank test and HR p values <0.05).
Eukaryotic translation initiation factors were introduced as novel drug targets in many cancers [59]. Eukaryotic translation initiation factor 2 subunit 1 (EIF2S1 [EIF2A]) mediates the binding of Met-tRNAi to the 40S/mRNA complex [60]. Several reports have linked the missregulation of eukaryotic translation initiation factors to abnormal cell growth and tumor development [61][62][63]. Additionally, angiogenesis and metastasis in cancer cells are linked to the PERK/eIF2a signaling pathway [64,65]. Intestinal-type adenocarcinoma is rare cancer afecting the nasal cavity and paranasal sinuses. Schatz et al. [66] found that it signifcantly overexpressed EIF2S1, EIF5A, and EIF6 when compared to healthy tissue samples. Recently, Li et al. [67] linked the UTP14A overexpression in oesophageal squamous cell carcinoma (ESCC) cells to the upregulation of PERK/eIF2a signaling pathway, leading to the cell cycle process and migration of ESCC cells. Te overexpression of    [68] performed a study to elucidate the role of caveolin-1 (CAV-1) and ferroptosis on HNSCC development. Te research fndings by Lu et al. [68] showed that CAV-1 was markedly overexpressed in HNSCC compared to healthy tissues. Te ferroptosis process was signifcantly inhibited by CAV-1, which increased cell proliferation, invasion, and metastasis. CAV-1 was signifcantly associated with a dismal prognosis in HNSCC patients as well. Terefore, CAV-1 was assigned as a potential target in HNSCC.
Te present study had several limitations. Te dataset GSE168227 only included patients from the Head and Neck Clinic of Regional Cancer Centre (Tiruvananthapuram, India). As a result, patients from other nations could not       completely beneft from the current results. Additionally, the sample size was small since the GSE168227 only included 48 oral tissue samples. A set of miRNAs selected in the dataset GSE168227 was based on the GPL8227 platform, which may not symbolize all the miRNAs.

Conclusion
Te present study suggests has-miR-136 and has-miR-377 as common DEMs in OTSCC and BSCC compared to the normal oral mucosa (p value <0.01; |Log2 FC| > 1). A total of 976 genes were identifed as validated targets of common DEMs. Ninety-six genes showed salient centrality within the PIM mediating the tumorigenesis of OTSCC and BSCC, in which EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, and HSPA5 were signifcantly linked to a poor prognosis in HNSCC. In addition, individuals with HNSCC had a better prognosis when they had NTRK2, HNRNPH1, DDX17, and WDR82 overexpression. "Clathrin-mediated endocytosis" was signifcantly enriched in OTSCC and BSCC. Our fndings might improve the prognosis of patients with OTSCC/BSCC, leading to more efective therapeutic strategies. However, in vitro and in vivo confrmations are needed in the future.

Data Availability
Te datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Ethical Approval
Te present study was approved by the Ethics Committee of the Hamadan University of Medical Sciences, Hamadan, Iran (ethics no. IR.UMSHA.REC.1401.451).

Conflicts of Interest
Te authors declare that there are no conficts of interest.