Abstract

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 ( 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.

1. 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 affected by lip and oral cavity cancer in 2020, leading to 177,757 deaths [1]. The 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 classified 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) [58]. 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 efficiently. 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 identified after they are already advanced [11]. In order to battle OSCC and increase patient survival rates, it is necessary to suggest more effective 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 [1315]. 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]. Therefore, miRNAs showed regulatory roles in critical biological procedures, including the cell cycle process, apoptosis and necroptosis, proliferation, and differentiation [18]. Consequently, miRNAs were defined 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 identified 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.

The present study hypothesized that miRNAs might participate in the tumorigenesis of OTSCC and BSCC. It was suggested that there are common differentially expressed miRNAs (DEMs) in the OTSCC and BSCC tissues compared to the normal oral mucosa, leading to abnormal expression of their target genes. These 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. Thus, utilizing the orthogonal partial least squares (OPLSs), common DEMs in OTSCC and BSCC in comparison to their respective healthy tissues were revealed. Then, validated targets of common DEMs were identified, 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 identified.

The 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. The 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 (Thiruvananthapuram, 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. The present study was confirmed by the Ethics Committee of Hamadan University of Medical Sciences, Hamadan, Iran (ethics no. IR.UMSHA.REC.1401.451).

2. Methods

2.1. Dataset Recovery and Statistical Analysis

The miRNA expression dataset GSE168227 [23], based on the platform GPL8227 (Agilent-019118 Human miRNA Microarray 2.0 G4470B), was downloaded as a TXT file from Gene Expression Omnibus (GEO), available at https://www.ncbi.nlm.nih.gov/geo [24]. The 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-off condition was set to 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.

2.2. PPI Network Analysis Based on Common DEM Targets

The experimentally validated targets of common DEMs were identified 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.

2.3. Gene Set Enrichment Analysis

Significant molecular pathways and Gene Ontology (GO) terms involved in the malignant transformation of normal oral mucosa to OTSCC and BSCC were unraveled using g:Profiler tool, available at https://biit.cs.ut.ee/gprofiler/gost [34]. The cut-off condition was set to the false discover rate (FDR) < 0.05 and the number of entities ≥10, following our previous study [35]. The 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].

2.4. 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 Profiling Interactive Analysis 2 (GEPIA2) database [38], available at http://gepia2.cancer-pku.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 defined as genes with the log-rank test and the hazard ratio (HR) values 0.05. The prognostic role of the combination of genes was also evaluated [39].

2.5. Gene Expression Evaluation of Prognostic Markers

The 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].

3. Results

3.1. DEMs in OTSCC and BSCC

Two prediction models were created using OPLSs to find DEMs between the OTSCC and BSCC tissues in comparison to their healthy counterparts. Each model was built using 369 variables and 24 samples. At 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 identified 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 value and Log 2 FC (Figure 1(b)).

3.2. Hub Genes, Modules, Pathways, and GO terms

The TarBase database detected 976 genes as experimentally validated targets for common DEMs. The 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. The average degree, betweenness, and closeness value was calculated as 14.33, 0.0023, and 0.33, respectively. Ninety-six nodes were then identified as hub genes linked to the pathogenesis of OTSCC and BSCC. (Table 2). Further structural analysis was performed using the MCODE plugin. Seven substantial clusters (cluster No. 1, cluster No. 2, cluster No. 3, cluster No. 4, cluster No. 7, cluster No. 8, and cluster No. 9) were found to be involved in the pathways and BPs linked to the pathogenesis of OTSCC and BSCC. Cluster No. 1 demonstrated the most MCODE score (MCODE score = 17.714), and cluster No. 2 included the most number of genes (n = 55) (Figure 2). The 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). The 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 affected 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 significant pathways and GO terms dysregulated in the pathogenesis of OTSCC and BSCC.

3.3. Prognostic Markers

Kaplan–Meier curves showed that the overexpression of EIF2S1, CAV1, RAN, ANXA5, CYCS, CFL1, MYC, HSP90AA1, PKM, and HSPA5 was significantly 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 value = 0.00016, and HR 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 value, and HR value of 0.71, 0.011, and 0.011, respectively, NTRK2 was the most positive marker. Therefore, 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. The combination of EIF2S1, CAV1, RAN, and ANXA5 showed a considerable negative signature with the criteria of HR = 1.6, log-rank value = 0.00028, and HR value = 0.00032 (Table 3).

3.4. Gene Expression of Markers

Based on boxplot analysis, the mRNA levels of EIF2S1, CAV1, RAN, SELE, ANXA5, CYCS, CFL1, HSP90AA1, PKM, HSPA5, HNRNPH1, and DDX17 exhibited a significant overexpression in HNSCC than in healthy controls. NTRK2 demonstrated a considerable underexpression in HNSCC than in normal oral mucosa. The MYC and WDR82 expressions were insignificant (Figure 6).

4. Discussion

The incidence of OTSCC and BSCC is considered the first and second highest among all oral cancers, respectively [40, 41]. MMP-13 inhibitors [42] may be effective 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) significantly 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 influences the pathways that lead to those receptors’ downstream signals [4345]. Xiao et al. [46] showed that ERK1/2 phosphorylates the FCH/F-BAR and SH3 domain-containing 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.

The 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 significant 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, ), suggesting the tumor suppressive role of miR-377-3P by downregulating the E2F3. The 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 significant 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 [5052]. Thiebaut et al. [53] showed a negative correlation between the expression of the miR-136-5p and ERα36 in breast cancer cells. The 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 effect 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 effect of the genes in patients with HNSCC. Fang et al. [57] identified 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. Then, Dai et al. [58] used the GEPIA for HNSCC to conduct a Kaplan–Meier survival analysis to investigate the possible predictive significance of the numerous hub CpG site-associated genes. The 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 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 [6163]. Additionally, angiogenesis and metastasis in cancer cells are linked to the PERK/eIF2a signaling pathway [64, 65]. Intestinal-type adenocarcinoma is rare cancer affecting the nasal cavity and paranasal sinuses. Schatz et al. [66] found that it significantly 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. The overexpression of EIF2S1 was confirmed at the mRNA level in HNSCC compared with healthy tissues.

Lu et al. [68] performed a study to elucidate the role of caveolin-1 (CAV-1) and ferroptosis on HNSCC development. The research findings by Lu et al. [68] showed that CAV-1 was markedly overexpressed in HNSCC compared to healthy tissues. The ferroptosis process was significantly inhibited by CAV-1, which increased cell proliferation, invasion, and metastasis. CAV-1 was significantly associated with a dismal prognosis in HNSCC patients as well. Therefore, CAV-1 was assigned as a potential target in HNSCC.

The present study had several limitations. The dataset GSE168227 only included patients from the Head and Neck Clinic of Regional Cancer Centre (Thiruvananthapuram, India). As a result, patients from other nations could not completely benefit 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.

5. Conclusion

The present study suggests has-miR-136 and has-miR-377 as common DEMs in OTSCC and BSCC compared to the normal oral mucosa ( value <0.01; |Log2 FC| > 1). A total of 976 genes were identified 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 significantly 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 significantly enriched in OTSCC and BSCC. Our findings might improve the prognosis of patients with OTSCC/BSCC, leading to more effective therapeutic strategies. However, in vitro and in vivo confirmations are needed in the future.

Data Availability

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

Ethical Approval

The 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

The authors declare that there are no conflicts of interest.

Authors’ Contributions

AT and S.Sh performed the study design. AT conducted the statistical analysis. The PIM, clustering, gene-set enrichment, and survival analyses were managed and illustrated by all authors. The results were analyzed and discussed by AT, S.Sh, and P.M. AT wrote the paper. S.Sh and Sh.J edited the paper. The final version of the article was read and approved by all authors.

Acknowledgments

The authors would like to thank the Research Center for Molecular Medicine and the Dental Research Center, the Hamadan University of Medical Sciences, Hamadan, Iran, for their support.