Analysis of the Differentially Expressed Genes Induced by Cisplatin Resistance in Oral Squamous Cell Carcinomas and Their Interaction

Background Oral squamous cell carcinoma (OSCC) is a solid tumor, which originates from squamous epithelium, with about 400,000 new-cases/year worldwidely. Presently, chemoradiotherapy is the most important adjuvant treatment for OSCC, mostly in advanced tumors. However, clinical resistance to chemotherapy still leads to poor prognosis of OSCC patients. Via high-throughput analysis of gene expression database of OSCC, we investigated the molecular mechanisms underlying cisplatin resistance in OSCC, analyzing the differentially expressed genes (DEGs) and their regulatory relationship, to clarify the molecular basis of OSCC chemotherapy resistance and provide a theoretical foundation for the treatment of patients with OSCC and individualized therapeutic targets accurately. Methods Datasets related to “OSCC” and “cisplatin resistance” (GSE111585 and GSE115119) were downloaded from the GEO database and analyzed by GEO2R. Venn diagram was used to obtain drug-resistance-related DEGs. Functional enrichment analysis and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis were performed on DEGs using The Database for Annotation, Visualization and Integrated Discovery (DAVID) software. Protein–protein interaction (PPI) network was constructed by STRING (search tool for recurring instances of neighbouring genes) database. Potential target genes of miRNA were predicted via miRDB, and cBioportal was used to analyze the function and survival of the potential functional genes. Results Forty-eight upregulated DEGs and 49 downregulated DEGs were obtained from the datasets, with cutoff as p < 0.01 and |log FC| > 1. The DEGs in OSCC mainly enriched in cell proliferation regulation, and chemokine activity. In PPI network with hub score > 300, the hub genes were identified as NOTCH1, JUN, CTNNB1, CEBPA, and ETS1. Among miRNA–mRNA targeting regulatory network, hsa-mir-200c-3p, hsa-mir-200b-3p, hsa-mir-429, and hsa-mir-139-5p were found to simultaneously regulate multiple hub genes. Survival analysis showed that patients with high CTNNB1 or low CEBPA expression had poor outcome. Conclusions In the OSCC cisplatin-resistant cell lines, NOTCH1, JUN, CTNNB1, CEBPA, and ETS1 were found as the hub genes involved in regulating the cisplatin resistance of OSCC. Members of the miR-200 family may reverse drug resistance of OSCC cells by regulating the hub genes, which can act as potential targets for the treatment of OSCC patients with cisplatin resistance.


INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC), the sixth most common malignant tumor in the world (Kim et al., 2011a), is an important public health issue worldwide. Among the total HNSCC cases, 30% are oral squamous cell carcinoma (OSCC) cases (World Health, 2003;Petersen, 2003a;Petersen, 2003b). In 2012, about 145,000 patients with OSCC died worldwide, with a mortality rate of 1.8% (Petersen, 2005;Kim et al., 2011a;Ong et al., 2016). Interestingly, OSCC is one of the three most common malignancies in Central and South Asia. In India, the age-standardized incidence of OSCC is 12.6 per 100,000 people (Petersen, 2005). According to statistics, the incidence of OSCC has increased sharply in several countries and regions, including Denmark, France, Germany, Scotland, and Central and Eastern Europe (Petersen, 2005).
OSCC can occur in different areas of the mouth and tongue, including lips, alveolar ridge, oral floor, oral tongue, hard palate, posterior molars triangle, and buccal mucosa, lined by squamous epithelium and scattered in smaller salivary glands and lymphatic drainage pathways. OSCC is common in the elder people with a history of tobacco and alcohol usage, with malignant tumors or somatic cell mutation by inducing DNA damage (Leemans et al., 2011). Although surgery is the main treatment strategy for OSCC, chemoradiotherapy is also an effective method, especially for advanced tumors. However, drug resistance due to unraveled molecular mechanisms significantly reduces the survival of OSCC patients.
Since the first miRNA-lin-4 was identified in 1993, miRNAs have attracted the attention of researchers in the field of gene expression regulation and gene therapy (Liang et al., 2014). By inhibition of RNA translation or degradation of target mRNA, miRNAs act as negative gene regulators at the posttranscriptional level (Sakai et al., 2013). Importantly, miRNAs can simultaneously modulate many target genes, such as tumor suppressors or oncogenes, widely influencing the phenotype of malignant tumors. Since miRNAs have been found to have important role in various aspects of malignant tumors, including oncogenesis, proliferation, metastasis, multidrug resistance, selfrenewal, and differentiation of malignant stem cells (Wu et al., 2014), they may represent a new set of therapeutic target biomarkers for finding multidrug resistance in malignant tumors (Hong et al., 2013).
In this study, the potential molecular mechanisms of cisplatin resistance of OSCC were studied by using high-throughput gene expression database. The differentially expressed genes (DEGs) in OSCC and their regulatory relationships were analyzed, in order to elucidate the molecular basis of OSCC chemotherapy resistance, and to provide theoretical basis and individualized precise therapeutic targets for the treatment of OSCC patients.

Microarray Datasets
"OSCC" and "cisplatin resistance" were used as the keywords for searching the GEO database, and GSE111585 and GSE115119 were downloaded as the gene expression data sets for cisplatin resistance in OSCC; the platforms used were GPL14715 and GPL16955.
GSE111585 included six samples of SCC9 cells and was divided into normal group and drug resistance group (Lin et al., 2018). GSE115119 contained four samples of CAL-27 cells, with normal group and drug resistant group. Both SCC9 and CAL-27 are human OSCC cell lines.

Data Analysis and Differential Expressed Gene Acquisition
Limma package of R software (GEO2R) was used for analysis of the original datasets. |log FC| > 1 and p value < 0.01 were defined as the cutoff values for further analysis of DEGs. Volcano maps were constructed by SangerBox software.

Functional Enrichment Analysis of DEGs
Gene Ontology (GO) provides a computational model of biological systems, from the molecular to the organism level, across different species in the following three categories: biological process (BP), molecular function (MF), and cellular component (CC) (Thomas, 2017). Kyoto encyclopedia of genes and genomes (KEGG) is a database for high-level functions and utilities of the biological systems, based on molecular-level information of genome sequencing and other high-throughput experimental technologies (Kanehisa et al., 2017). DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/) comprises a comprehensive set of functional annotation tools for functional enrichment analysis of gene groups (Huang da et al., 2009). To identify the biological significance of DEGs in cisplatin-resistant OSCC cells, DAVID 6.8 was used to analyze GO function and KEGG pathway enrichment, with the enrichment standard as p < 0.05.

Protein-Protein Interaction Network of DEGs
Protein-protein interaction (PPI) network analysis is helpful to investigate the molecular mechanisms of diseases and discover new drug targets from a systematic perspective. STRING 11.0 (https://string-db.org/), covering more than 5,000 organisms with known and predicted protein-protein interactions, provides direct (physical) and indirect (functional) association (Szklarczyk et al., 2019). The PPI analysis of DEGs was performed by STRING 11.0, and the results were analyzed by Cytoscape 3.7.1. Furthermore, the cytoHubba plug in was used to calculate the interaction coefficient score between the DEGs. The top genes with hub score > 300 were identified as the hub genes with high connectivity in the PPI network.

Predicting Hub Gene-Related miRNAs
MicroRNAs (miRNAs), small non-coding RNA molecules with highly conserved regions, regulate the expression of target genes by binding to the 3'-untranslated regions (3'-UTR) of specific mRNAs, involved in many physiological and disease processes. Each miRNA is thought to regulate multiple genes with enormous potential regulatory circuitry afforded by miRNA (Lim et al., 2003). To identify the potential miRNA-mRNA interaction in the network of the hub genes, miRDB (http:// mirdb.org/), an online resource for miRNA target prediction and functional annotation (Wong and Wang, 2015), was used to predict the hub gene-related miRNAs, and the miRNA-mRNA regulatory network was constructed by Cytoscape 3.7.1.

Expression and Survival Analysis of Hub Genes
The Oncomine database (https://www.oncomine.org/resource/ login.html), an online cancer microarray database-mining platform (Rhodes et al., 2004), was used to investigate the difference in transcriptional levels of the hub genes in HNSCC and normal tissues.
As mutations of oncogenes and/or tumor-suppressor genes are frequent in tumor tissues, the Human Protein Atlas (http://www.proteinatlas.org/) was analyzed for the prognostic values of the hub genes (Uhlen et al., 2017), and cBioportal database (http://www.cbioportal.org/), an open-access online resource for multi-dimension analysis of data from The Cancer Genome Atlas (TCGA) (Gao et al., 2013), was used to analyze the effects of mutations in hub genes on the survival of patients with OSCC (MD Anderson, Cancer Discov, 2013).

Difference of Gene Expression Between Parental and Cisplatin-Resistant OSCC Cells
The gene expression microarray datasets, GSE111585 and GSE115119 were downloaded from GEO datasets with paired parental and cisplatin-resistant OSCC cells. As shown in Figure 1, the expression of most genes in cisplatin-resistant OSCC cells was similar to that of the parental OSCC cells. Cluster analysis by R software (|log FC| > 1 and p value < 0.01 as the cutoff) revealed 1,386 up-regulated genes and 643 down-regulated genes in cisplatin-resistant OSCC cells compared with parental OSCC cells in GSE111585 ( Figure 1A), and 757 up-regulated genes and 625 down-regulated genes in cisplatin-resistant OSCC cells compared with parental OSCC cells in GSE115119 ( Figure 1B).
The intersection between DEGs and the list of oncogenes drawn by Venn software showed 48 up-regulated DEGs ( Figure 1C), and 49 down-regulated DEGs were obtained via intersection between down-regulated genes and the list of tumorsuppressor genes ( Figure 1D).

Close Association of the DEGs With the Regulation of Transcription and microRNAs in Cancers
Using the DAVID analysis software, functional enrichment analyses (BP, MF, and CC) of the DEGs were done. BP enrichment showed that the up-regulated DEGs were mainly enriched in cell proliferation regulation, inflammatory reaction, lipopolysaccharide, cells in response to growth factors to stimulate, neuronal migration, transmembrane receptor protein tyrosine kinase signaling pathway, and transcription of RNA polymerase II promoter (Figure 2A), whereas down-regulated DEGs were significantly enriched mainly in the following GO terms: response to X-ray, RNA polymerase II promoter negative transcription regulation, and DNA damage response ( Figure 2B).
MF enrichment showed that the up-regulated DEGs were significantly enriched in chemokine activity, transcription factor activity, sequence specific DNA binding, non-membrane crossing protein tyrosine kinase activity, and sequence specific DNA binding ( Figure 2C), and the down-regulated DEGs were enriched in p53 binding, sequence specific DNA binding, transcriptional activator activity, and RNA polymerase II hub promoter proximal region sequence specific binding ( Figure 2D).
CC analysis predicted close association between the upregulated DEGs and the following GO terms: mRNA cutting, polyadenylation specific factor complex, extracellular space, promyelocytic leukemia proteome, and transcription factor complex ( Figure 2E), and significant relation was found between the down-regulated DEGs and the following GO terms: junction complex, desmosomes, ciliated tips, cytoplasm, nuclear cytoplasm, and plasma membrane ( Figure 2F).
KEGG pathway analysis provided the potential function cluster of DEGs, showing that the up-regulated DEGs were clustered in malaria, human T-cell leukemia virus type I, the way of malignant tumor, legionella infection disease, TNF signaling pathways, and T-cell receptors signaling pathways ( Figure 3A), whereas the down-regulated DEGs were significantly concentrated in axon guidance and microRNAs in cancers ( Figure 3B).

Identification of Hub Genes Through PPI Network of DEGs
To further analyze the correlation between DEGs in cisplatinresistant OSCC cells, STRING was used to construct PPI network showing close relationship between the DEGs (Supplemental Figure 1), and their hub score was calculated. The genes with high hub score were predicted to have a strong association with other genes (shown in dark color in the figures). As shown in Figure 4, based on the cutoff hub score > 300, the following five genes were selected as the hub genes: NOTCH1, JUN, CTNNB1, CEBPA, and ETS1.

Construction of miRNA-mRNA Network Based on Predicting miRNA-Target Genes
As the DEGs in cisplatin-resistant OSCC cells were closely related to tumor-related miRNA, miRDB database was used to predict potential miRNAs that might participate in the transcriptional regulation of the hub genes in this process. The prediction scores were also collected from the miRDB database, and the miRNA-mRNA with high score meant close potential function of miRNA in regulation of the target mRNA. After setting cutoff > 80, Cytoscape software was used to construct the miRNA-mRNA network ( Figure 5). Interestingly, hsa-miR-200c-3p, hsa-miR-200b-3p, hsa-miR-429, and hsa-miR-139-5p could simultaneously regulate multiple hub genes, which may be the key miRNAs involved in this process. Interestingly, hsa-miR-200c-3p, hsa-miR-200b-3p, and hsa-miR-429 belong to miR-200 family members, with similar functions; suppression of ZEB1/2, followed by inhibition of epithelial-mesenchymal transition (EMT).

The Expression Pattern of Hub Genes in OSCC
To investigate the potential function of the hub genes in OSCC, Oncomine database was used to analyze the difference in the expression levels of the hub genes. However, due to limited research on OSCC, only one study revealed that the expression of CTNNB1 and ETS1 in tumor tissues was higher than that in normal tissues, with 2.285 and 2.111 fold change, respectively, while no difference was found in the expression of NOTCH1 and JUN genes in the two tissues ( Figure 6).

Survival Value of Hub Genes in OSCC
For survival analysis, cBioportal based on TCGA database was used, which revealed that low expression of CTNNB1 in patients with OSCC showed better overall survival (p = 0.01) ( Figure 7C), and low expression of CEBPA predicted poor overall survival in  OSCC patients (p = 0.04) ( Figure 7D). Although the expression of other hub genes did not show a significant relationship with the survival status of OSCC patients (p > 0.05), the OSCC patients with high expression of NOTCH1 ( Figure 7A) and ETS1 ( Figure 7E) or low expression of JUN ( Figure 7B) tended to have long lifespans. The dataset obtained from MD Anderson, Cancer Discov 2013, showed that the median overall survival of all OSCC patients was 78.8 months. Except for NOTCH1, no mutation was found in the other hub genes in the OSCC patients. And the mutations in NOTCH1 showed no significant association with the overall survival of patients with OSCC (Figure 8), suggesting that the regulation, without mutation of the hub genes was the main mechanism of cisplatin resistance in OSCC.

DISCUSSION
Worldwide, OSCC is an important public health issue with limited therapy strategies and researches; systemic drug resistance has aggravated this situation. In this study, highthroughput screening was used to explore the potential genes involved in cisplatin resistance of OSCC, and NOTCH1, JUN, CTNNB1, CEBPA, and ETS1 were identified as the hub genes in the occurrence of cisplatin resistance. These genes were found to be regulated by the members of the miR-200 family. Regulation of the corresponding hub genes by miRNAs may reverse cisplatin resistance of OSCC, and the sensitivity of tumor cells to cisplatin maybe restored; thus, providing a novel potential target for anticancer therapy.
Studies have shown that changes in NOTCH signaling pathway are associated with many human cancers (Villanueva et al., 2012). NOTCH1 is reported to be both a tumor suppressor gene and a tumor oncogene. The tumorigenic or anti-tumor activity of NOTCH family members in different types of tumors displays its role in promoting or inhibiting the undifferentiated state of stem cells in the corresponding tissues (Wang et al., 2012). Carcinogenic action of NOTCH has been found in many cancers, including non-small cell lung cancer (Lenhart et al., 2015), acute T lymphoblastic leukemia (Weng et al., 2004), and malignant gliomas (Purow et al., 2005). In contrast, NOTCH1 signaling is inhibited in neuroendocrine tumor cells, including small cell lung cancer (Platta et al., 2008). This suggests that induction of NOTCH1 expression is an effective strategy for treating these tumors. NOTCH signaling pathway is also involved in chemotherapy resistance. For example, NOTCH1 plays an important role in cisplatin resistance mechanism of head and neck squamous cell tumor, colorectal tumor, ovarian cancer (Wang et al., 2010), and other malignant tumors. In this study also, expression of NOTCH1 gene was found to be significantly inhibited in cisplatin-resistant OSCC cell lines as compared to that in normal or tumor tissues, but no effect was observed on the overall survival of patients. These results suggest that NOTCH1 signaling molecules may be involved in different FIGURE 4 | The PPI network of DEGs with Hub score. The dark color indicates high hub score, and the light color predicts low hub score.
biological processes of malignant tumor development through different molecular pathways, and could play an important role in resistance of OSCC against cisplatin and other chemotherapy drugs.
JUN is a protein-coding gene and has no introns; it is located in 1p32-p31: a chromosomal region involved in human malignant translocations and deletions (Fazal et al., 2017) JUN-related diseases, include sarcomas and whooping cough (Syc-Mazurek et al., 2017). JUN is involved in the following pathways: apoptosis regulation, signal transduction, tacrolimus/ cyclosporine pathway, and pharmacodynamics. JUN is also associated with sequence-specific DNA binding (GO annotation). In this study, although JUN molecular expression was significantly changed, its correlation with malignant tumor tissues and its influence on patient survival were not found. Therefore, its function and molecular mechanism will be explored in future studies.
The protein encoded by CTNNB1 is a part of the protein complex that forms the adhesive junctions. Adhesion is necessary to create and maintain the epithelial layer (Li et al., 2017a). The coding proteins, which also include the actin cytoskeleton, are responsible for signaling contact inhibition, and once the upper cortex completes signaling, the cell stops dividing. Finally, the   protein binds to the product of the APC gene, which is mutated in colorectal adenomatous polyposis. The mutation is a cause of colorectal cancer, hairy tumors, medulloblastoma, and ovarian cancer (Li et al., 2017b). Selective splicing of CTNNB1 RNA leads to multiple transcript variants. Diseases associated with CTNNB1, include hairy tumors and intellectual disability, both being 19 autosomal dominant (Lee et al., 2018). The pathways associated with CTNNB1 are beta-adrenergic signaling and blood-brain barrier pathways. Because it inhibits the expression of downstream signals, GO annotations associated with it, include DNA-binding transcription factor activity and binding. In this study, it was found that the expression of CTNNB1 in tumor tissues was significantly higher than that in normal tissues, and the survival period of patients with high expression of CTNNB1 was significantly shortened. These results suggest that CTNNB1 also plays an important role in the occurrence and development of OSCC, but the mechanism of its influence on cisplatin chemotherapy resistance needs to be further studied and explored. CEBPA is an intron-free transcription factor that contains a basic leucine zipper domain and recognizes the CCAAT motif in the target gene promoter (Mannelli et al., 2017). The coding proteins act in homodimers and heterodimers with CCAAT or enhancer binding proteins, beta and gamma. The activity of CEBPA protein regulates the expression of genes, which are involved in cell cycle regulation and weight balance. Mutation in the CEBPA gene has been linked to acute myeloid leukemia (Avellino et al., 2016). CEBPA mutations are particularly associated with cytogenetically-normal AML (Taskesen et al., 2011). CEBPA is necessary for granulocyte formation in mice. Mutations in CEBPA are associated with longer survival of OSCC patients. CEBPA-related diseases, include leukemia, acute myeloid leukemia, and myeloid leukemia. The pathways associated with CEBPA are adenoid cystic carcinoma and the adipogenesis pathway. CEBPA has important DNA binding transcription factor activity and can bind to sequence specific DNA. However, there are no relevant studies on CEBPA and cisplatin resistance of OSCC at present. In the current study, we found that high expression of CEBPA is closely related to the prognosis of OSCC patients.
ETS1 is a member of the encoding transcription factor ETS family, which has a conserved DNA binding domain of ETS that recognizes the hub consistent DNA sequence GGAA/T in the target gene (Poon and Kim, 2017). These proteins act as transcriptional activators or inhibitors of many genes and are involved in stem cell development, cell aging and death, and tumorigenesis. Splicing transcriptional variants encoding different subtypes have also been previously described. Jacobsen syndrome and estrogen receptor negative breast cancer are the diseases associated with ETS1 (Carpinelli et al., 2015). The pathways involving ETS1 include focal adhesion and focal adhesion kinase mediated signal transduction events. The gene also has important DNA-binding transcription factor activity and transcription factor binding. We found that ETS1 is an important cisplatin resistant gene based on highthroughput data analysis, PPI network, and expression verification. Studies have shown that overexpression of ETS1 induces IKK alpha mRNA and protein expression as well as IKK alpha activity (Gu et al., 2004). In a previous study, ETS1 protein expression and IKK alpha were significantly upregulated in 231 cisplatin-resistant cell lines. Inhibition of ETS1 expression has been reported to enhance cisplatin sensitivity of resistant cell lines. ETS1 knockout increases the stability of cisplatin in mouse xenograft models . These results are similar to the results obtained in the current study. ETS1 was highly expressed in cisplatin-resistant OSCC cell lines as compared to that in the normal tissues; ETS1 was highly expressed in tumor tissues, suggesting that it is an important molecule in this process.
Based on previous studies on hub genes and members of the miR-200 family, miR-200b/a/429 transcription is known to be regulated by different transcriptional factors in tissue-specific manner (Kim et al., 2011b). ZEB1/2 is the classical target gene of miR-200s, and many other potential factors have also been reported as the genes regulated by miR-200s (Nagalla et al., 2011). In the current study, new potential target genes were reported as the hub genes in cisplatin-resistant OSCC cells. In 2018, Liu et al. reported a smart miRNA-reporter gene for in vitro and in vivo imaging of biogenesis of miRNA and their related functions (Liu et al., 2018). Further study involving the reporter system could be helpful in investigation of the relationship between miR-200s and the hub genes in OSCC.
And as the researches related to OSCC are limited, the relationship between the expression of hub genes and clinicopathological parameters in OSCC patients will be collected and analyzed in the further, to confirm their roles in the occurrence of cisplatin resistance in OSCC.

CONCLUSION
We found that NOTCH1, JUN, CTNNB1, CEBPA, and ETS1 were the key genes regulating cisplatin resistance in OSCC drugresistant cell lines, and the miR-200 family may be capable of reversing OSCC cell resistance by regulating NOTCH1, JUN, and ETS1, which could also act as potential targets for treating cisplatin resistant OSCC patients.

DATA AVAILABILITY STATEMENT
GSE111585 and GSE115119 were downloaded from Gene Expression Omnibus.

AUTHOR CONTRIBUTIONS
JL conceptualized and designed the study. W-TC, H-TW, and JL organized the database, searched literature, and structured and drafted the manuscript. G-WL, J-XS, Q-QY, M-LZ, and W-JC analyzed and interpreted the data, and drafted and revised the manuscript. JL revised the original manuscript. All authors contributed to manuscript revision and read and approved the submitted version.