Integrated bioinformatics analysis of potential biomarkers and candidate drugs of esophageal squamous cell carcinoma

Ming-Qian Hao ,Ting Zhao ,Da-Yong Chen ,Feng-Xian Zhang ,Yi-Wen Zhang ,Qi-Teng Ding ,Shu-Wen Sun,Jin-Ping Zhang,Ling Dong,Chuan-Bo Ding,Wen-Cong Liu*

1College of Traditional Chinese Medicine,Jilin Agricultural Science and Technology College,Jilin 132109,China.2School of Chinese Medicinal Materials,Jilin Agricultural University,Changchun 130118,China.3Jilin Province Center for Drug Inspection,Changchun 130062,China.

Abstract Esophageal squamous cell carcinoma (ESCC),the major subtype of esophageal carcinoma(ESCA),is one of the most lethal malignancies worldwide.This study aimed to identify potential biomarkers and/or therapeutic targets for ESCC.The datasets GSE44021,GSE77861,GSE20347,and GSE29001 retrieved from the Gene Expression Omnibus (GEO)database contained 117 ESCC tissues and 109 normal tissues.Differentially expressed genes(DEGs) associated with ESCC were identified using the GEO2R tool.Dysregulated pathways associated with ESCC mainly included mitotic regulation,cell cycle,ECM-receptor interaction,DNA replication,etc.The protein-protein interaction (PPI) network of overlapping DEGs was constructed and nine key genes (KGs) were identified from the complex interaction network using Degree,maximum neighborhood component (MNC),and maximal clique centrality(MCC)algorithms.Expression patterns of KGs at the transcriptional and translational levels were validated using ESCC-related data from the Cancer Genome Atlas (TCGA),Oncomine,and Human Protein Atlas (HPA) databases.Genetic alterations calculation,immune cell infiltrates evaluation,methylation analysis,prognostic analysis,transcription factors (TFs) and miRNAs regulatory networks construction,and targeted drug prediction were further performed.It was also found that the knockout of these KGs affected the survival of more than two types of ESCC cells by genome-wide CRISPR-Cas9 dropout screens.In conclusion,we identified KGs,TFs,and miRNAs with biomarker potential (e.g.,NDC80,BUB1,TOP2A,AURKA,AURKB,TTK,UBE2C,TPX2,BUB1B,E2F1,and hsa-miR-483-5p) and 23 candidate targeted drugs for ESCC by utilizing an integrated multi-omics approach.These findings provide additional insights into uncovering the molecular mechanism and improving the efficiency of clinical diagnosis and treatment for ESCC.

Keywords: esophageal squamous cell carcinoma;biomarker;transcription factor;miRNA;CRISPR-Cas9;candidate drug

Background

ESCC,the major pathological type of ESCA,is a global malignant tumor with poor diagnosis,prognosis,and survival rate [1].For the lack of effective preventive and diagnostic methods,most patients with ESCC are only diagnosed and detected at an advanced stage with a 5-year overall survival rate of less than 20% [2].In addition,although clinical therapies such as immunotherapy,neoadjuvant chemoradiotherapy,and minimally invasive esophagectomy have achieved some promising results,the risk of mortality from ESCC is still high due to population increase and aging,as well as incomplete treatment and tumor metastasis [3].Therefore,it is necessary to explore promising biomarkers and valuable therapeutic drugs for ESCC.

High-throughput sequencing technology and bioinformatics method have become powerful tools to monitor genome-wide gene expression changes and provide an indispensable platform for integrating public databases and exploring the molecular mechanisms of various human diseases [4,5].Wang et al.re-utilized the published transcriptome data to identify DEGs and signaling pathways associated with ESCC,and further found that aberrant expression of MAPK1,ACOX1,and SCP2 affected the overall survival time of ESCC patients[6].Han et al.demonstrated that AURKA was an abnormally methylated DEGs in ESCC,and its overexpression was associated with poor overall survival through an integrated bioinformatics analysis [7].Moreover,a recent multi-omics report identified 10 hub nodes of ESCC by expression regulatory network analysis of abnormal genes and investigated their potential as diagnostic and therapeutic targets for ESCC [8].Another bioinformatics research based on large-scale sequencing data reported that miR-139-3p was significantly down-regulated in the serum of patients with ESCC/ESCA and might be a candidate biomarker for predicting ESCC [9].The above bioinformatics research mainly focuses on whether single or multiple genes affect tumorigenesis and prognosis.However,the related work on gene function and stability,tumor immune microenvironment,and drug candidates is insufficient and needs to be further implemented.

Based on the above characteristics,we performed a comprehensive multi-omics task,including DEGs identification,PPI network construction,functional enrichment analysis,KGs identification,genetic alteration and methylation analysis,fitness gene prediction,tumor microenvironment and immune cell infiltration assessment,miRNAs-and TFs-KGs regulatory networks mapping,prognostic analysis,and targeted drug discovery,to reprocess large-scale sequencing data of ESCC from multiple human tumor-related databases.We believed this study could provide insights into the pathogenesis,clinical diagnosis,prognosis,and treatment of ESCC based on the identified biological mechanisms,biomarker targets,and drug candidates.

Methods

DEGs identification

Four gene expression profiles of ESCC (GSE44021,GSE77861,GSE29001,and GSE20347) were randomly searched from the GEO(https://www.ncbi.nlm.nih.gov/geo/) database and DEGs between ESCC tissues and normal tissues were identified by using GEO2R tool(https://www.ncbi.nlm.nih.gov/geo/geo2r/)with a screening criteria of adjustedP-value (adj.P.Val) ＜0.05 and |log2 fold change (logFC)|＞1 [10].Overlapping DEGs from different datasets were obtained using a VennDiagram package in R language (version 3.6.2) [11].Table S1 listed the details of the above datasets.

PPI network construction and KGs identification

PPI network was integrated using the STRING (version 11.0,https://string-db.org/) database with an interaction score of 0.7,and then visualized by the cytoscape software (version 3.7.2) [12-13].KGs were identified from the co-expression network using Degree,MNC and MCC algorithms of the CytoHubba plugin [14].

Expression verification,methylation assessment and prognostic analysis of KGs

The UALCAN (http://ualcan.path.uab.edu/index.html) and Oncomine(https://www.oncomine.org) databases were used to verify the transcriptional expression of KGs [15,16].Differences of these genes between normal and tumor tissues at the translational level were evaluated by the HPA (https://www.proteinatlas.org) database [17].We also assessed the methylation degree and prognostic value of KGs using the UALCAN tool.P＜ 0.05 was considered statistically significant.

Identification of fitness genes in ESCC cell lines

Fitness genes were defined as genes required for cell survival and growth.The Project Score database(https://score.depmap.sanger.ac.uk/),a public resource that contains genome-wide CRISPR-Cas9 dropout screening data in 323 cancer cell lines,was applied to evaluate fitness genes that promote the growth of three or more ESCC cell lines based on fitness scores [18].Corrected log fold change ＜-1 means that gene depletion prevents cell growth or viability.

Functional enrichment analysis

To elucidate the biological characteristics of overlapping DEGs,we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the Cluster Profiler package of R [19].Statistical significance was set atP＜0.05.

Evaluation of immune infiltration

Tumor Immune Estimation Resource (TIMER,https://cistrome.shinyapps.io/timer/) is an integrated platform for analyzing 32 cancer types’immune infiltration data[20].We explored the link between abnormal expression of KGs and the number of tumor-infiltrating immune cells (B Cell,CD4+T Cell,CD8+T Cell,Macrophage,Neutrophil,and Dendritic Cell) in ESCA samples.Correlation analysis between KGs was performed by the “Correlation”module.

Analysis of genetic alterations of KGs

cBioPortal (http://cbioportal.org) facilitates the exploration and visualization of complex cancer genomics data[21].In this study,227 samples from 2 ESCC reports in the cBioPortal database were used to explore genetic alterations in KGs.

Prediction of TFs and miRNAs-KGs interactions

To understand the transcriptional regulation process of KGs,we constructed a TFs-KGs interaction network using TRRUST platform(version 2,https://www.grnpedia.org/trrust/) [22].In addition,the miRNAs-mRNAs network was constructed to further explore the function of KGs and the regulation of miRNAs in ESCC.Briefly,we identified overlapping differentially expressed miRNAs (DEmiRNAs)from datasets GSE114110 and GSE43732 utilizing the GEO2R and VennDiagram tools with criteria of|logFC|＞1 and adj.P.Val ＜0.05.The interactions between KGs and overlapping DEmiRNAs were obtained from the miRDB (http://mirdb.org/),miRWalk(http://mirwalk.umm.uni-heidelberg.de/),and miRTarBase (http://mirtarbase.cuhk.edu.cn/php/index.php) databases [23-25].

Targeted drugs prediction

The DGIdb tool (http://www.dgidb.org/) can redesign and mine clinically relevant drug-gene interactions for personalized medicine[26].Hence,we further re-analyzed the compounds in the DGIdb database to identify anti-neoplastic drugs targeting KGs in ESCC.

Results

Identification of DEGs in ESCC

In the section on DEGs identification,we randomly selected four datasets involving a total of 117 ESCC samples and 109 paired normal samples.The data boxplot showed that the distribution of samples in each individual dataset tended to be consistent and could be used for further analysis (Figure S1).As presented in Figure 1A-D,1048 DEGs in GSE44021 (537 up-and 511 down-regulated),1260 DEGs in GSE77861 (684 up-and 576 down-regulated),4858 DEGs in GSE29001(2952 up-and 1906 down-regulated),and 1428 DEGs(660 up-and 768 down-regulated) in GSE20347 were identified.To avoid the bias of individual studies,we performed a VennDiagram analysis and found 462 overlapping DEGs (212 up-and 250 down-regulated)from the above datasets (Figure 1E and F,Table S2).

PPI network integration and KGs identification

A PPI network composed of overlapping DEGs was shown in Figure 2A,containing 255 nodes and 1151 interactions.Subsequently,nine KGs,namely NDC80,BUB1,TOP2A,AURKA,AURKB,TTK,UBE2C,TPX2,and BUB1B,were identified based on the Degree,MNC,and MCC algorithms (Figure 2B-E).Notably,these KGs were all significantly up-regulated in malignant lesion tissues and might be inseparable from the pathogenesis and progression of ESCC.

Expression verification,methylation assessment and prognostic analysis of KGs

Based on clinical data and pan-cancer reports retrieved from the TCGA and oncomine databases,we found that these KGs in ESCC were abnormally up-regulated in most human cancers (e.g.,Brain and CNS Cancer,ESCA,Breast Cancer,Lung Cancer,Bladder Cancer,and Sarcoma,Figure 2F and Figure S2A).

Figure 1 Identification of overlapping differentially expressed genes (DEGs) from four esophageal squamous cell carcinoma (ESCC)datasets.(A-D) Volcano plots of DEGs in GSE44021,GSE77861,GSE29001,and GSE20347,respectively.Red means up-regulated DEGs,green means down-regulated DEGs,and black means no significant difference.(E) Venn diagram of up-regulated DEGs.(F) Venn diagram of down-regulated DEGs.DEGs,differentially expressed genes;ESCC,esophageal squamous cell carcinoma.

Figure 2 Protein-protein interaction (PPI) network construction and key genes (KGs) identification.(A) PPI network of overlapping DEGs.Red diamonds: up-regulated genes;Blue ellipses: down-regulated genes.(B-D) Top 10 genes identified from the PPI network via MCC,MNC,and Degree algorithms,respectively.(E) Venn diagram for identifying KGs.(F) Heatmap of KGs expression in esophageal carcinoma (ESCA) data from the cancer genome atlas(TCGA)database.PPI,protein-protein interaction;KGs,key genes;ESCA,esophageal carcinoma;TCGA,the cancer genome atlas.

Further,we verified the expression trends of nine KGs in ESCC tissues by analyzing relevant samples from the TCGA database (Figure S2B-J).Interestingly,these genes were also over-expressed in esophageal adenocarcinoma (P＜0.001).We then compared the expression differences of KGs at the translational level in ESCC and normal tissues using immunohistochemical staining data from the HPA database.TOP2A is over-expressed in ESCC tissues and moderately expressed in normal tissues (Figure 3A and B).AURKA showed medium expression in tumor samples but was negative in adjacent normal tissues (Figure 3C and D).AURKB displayed medium expression in both ESCC and adjacent tissues (Figure 3E and F).Medium expression of TPX2 was found in the nucleus of ESCC cells and the cytoplasm/membrane of normal tissue cells (Figure 3G and H).UBE2C was found to have high expression in normal tissues but moderately expressed in ESCC tissues,which might be relevant to the relatively small sample size of ESCC included in the HPA database(Figure 3I and J).Unfortunately,we did not find the expression data of NDC80,BUB1,TTK,and BUB1B proteins in ESCC from the HPA database.In summary,we predict that the abnormal expression of these proteins is highly likely to participate in the pathogenesis and progression of ESCC,which needs further experimental confirmation.

Furthermore,we analyzed the methylation data of ESCC included in the TCGA database and found that the methylation levels of BUB1,UBE2C,and BUB1B were significantly higher in ESCC tissues than in normal tissues (P＜0.01 orP＜0.001);while the methylation levels of NDC80,BUB1,TOP2A,and TTK in esophageal adenocarcinoma were higher than those in normal tissues (P＜0.05 orP＜0.001,Figure S3).These findings indicated that KGs might be associated with the prognosis of ESCC,but current data showed that there was no significant relationship between the abnormal expression of these KGs and the overall survival of patients with ESCC (P＞0.05,Figure S4).Therefore,we recognized that further evidence was needed for prognostic analysis of ESCC.

Figure 3 Expression validation of key genes (KG) at the translational level by the Human Protein Atlas database.(A) TOP2A protein in normal tissues (staining: medium;intensity: strong;quantity: ＜25%;location: nuclear).(B) TOP2A protein in tumor tissues (staining: high;intensity:strong;quantity: 75%-25%;location:nuclear).(C)AURKA protein in normal tissues(staining:not detected;intensity:negative;quantity:none;location: none).(D) AURKA protein in tumor tissues (staining: medium;intensity: moderate;quantity: 75%-25%;location:cytoplasmic/membranous,nuclear).(E)AURKB protein in normal tissues(staining:medium;intensity:strong;quantity:＜25%;location:nuclear).(F) AURKB protein in tumor tissues (staining: medium;intensity: strong;quantity: ＜25%;location: nuclear).(G) TPX2 protein in normal tissues(staining: medium;intensity: moderate;quantity: 75%-25%;location: cytoplasmic/membranous).(H) TPX2 protein in tumor tissues (staining:medium;intensity: strong;quantity: ＜25%;location: nuclear).(I) UBE2C protein in normal tissues (staining: high;intensity: strong;quantity: ＞75%;location: nuclear).(J) UBE2C protein in tumor tissues (staining: medium;intensity: moderate;quantity: ＞75%;location: nuclear).KG,key genes.

Evaluation of fitness genes

The impacts of knockout KGs on the growth and cell viability of ESCC cell lines were further investigated using genome-scale CRISPR-Cas9 dropout screening data.In the results,all nine KGs could be defined as fitness genes of ESCC,that is,the deletion of these genes hindered the survival and growth of more than two ESCC cell lines.Knockout of TOP2A or AURKB,but not other KGs,most significantly inhibited the growth of 19 ESCC cell lines.In addition,KYSE-450,TE-10,and T-T cells,but not other ESCC cell lines,were more sensitive to KGs deletion,although there was no statistical significance (Figure 4).

GO and KEGG enrichment analysis

GO and KEGG enrichment analysis was further performed to elaborate the biological annotations of 462 overlapping DEGs in ESCC.GO enrichment mainly included three categories: biological process (BP),cell component (CC),and molecular function (MF) (Figure 5A).For the BP category,overlapping DEGs were mainly enriched in cell cycle-related terms,such as extracellular matrix organization,cell cycle checkpoint,mitotic nuclear division,and positive regulation of cell cycle,which was consistent with the biological features of the rapid proliferation of ESCC cells.CC enrichment analysis indicated that overlapping DEGs were related to the collagen-containing extracellular matrix,chromosomal region,endoplasmic reticulum lumen,and extracellular matrix component.In the MF category,overlapping DEGs were markedly involved in growth factor binding,extracellular matrix structural constituent,cell adhesion molecule binding,DNA replication origin binding,and so on.KEGG pathway analysis revealed that 462 overlapping DEGs were primarily linked to DNA replication,cell cycle,transcriptional mis-regulation in cancer,ECM-receptor interaction,protein digestion and absorption,and proteoglycans in cancer (Figure 5B).Furthermore,KGs such as BUB1B,TTK,and BUB1 mainly participated in the cell cycle.

These results were consistent with the understanding that abnormalities in cell growth regulatory factors and cell cycle were the main reasons for tumor formation,which further indicated that KGs and overlapping DEGs played crucial roles in the pathogenesis and tumor microenvironment of ESCC [27].

Characterization of immune infiltration and genetic changes of key genes

Tumor tissues contain not only tumor cells but also tumor-infiltrating immune cells closely related to tumor progression.In this study,nine KGs were correlated with purity and six types of immune-infiltrating cells (B Cell,CD4+T Cell,CD8+T Cell,Macrophage,Neutrophil,and Dendritic Cell) to varying degrees(Figure S5).Furthermore,there was a positive correlation in expression patterns among nine KGs in ESCC tissues,showing that these genes may be essential functional partners (Figure S6).

Figure 4 Effects of knockout of key genes (KGs) on cells of esophageal squamous cell carcinoma (ESCC).Corrected log fold change ＜-1 means that gene depletion prevented cell growth or viability.KGs,key genes;ESCC,esophageal squamous cell carcinoma.

Figure 5 Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis for overlapping differentially expressed genes (DEGs).(A) GO enrichment analysis of overlapping DEGs (Top 10 of each category).BP,biological process;CC,cellular component;MF,molecular function.The color represents P-value and the size of the spots represents the gene number.(B) KEGG pathway enrichment analysis.GO,gene ontology;KEGG,kyoto encyclopedia of genes and genomes;DEGs,differentially expressed genes.

Genetic alterations of 9 KGs in 227 ESCC patients were analyzed using the cBioPortal database.As a result,we observed that AURKA and BUB1B produced 0.9% alterations due to missense mutations in ESCC samples,respectively,and the genetic alteration rate of the remaining genes was 0% (Figure 6A).Therefore,these KGs were considered low-frequency mutant genes.

TFs-and miRNAs-KGs regulatory networks construction

TFs-KGs network containing 7 KGs (TOP2A,NDC80,AURKA,AURKB,TTK,UBE2C,and BUB1B) and 12 TFs(E2F1,E2F3,E2F4,MED1,etc.)was constructed (Figure 6B).We found that other TFs excluding ATF4 and OTX2 were significantly up-regulated in ESCC samples compared to normal tissues (P＜0.001,Figure S7).Interestingly,E2F1 not only activated AURKA and TOP2A genes but also affected the transcriptional process of AURKB.E2F3,MED1,and OTX2 as TFs all had potential activating effects on AURKA.In addition,ATF4 and CREB1 might activate the transcriptional process of NDC80,and ZNF143 and MYC had an activating effect on the expression of BUB1B and UBE2C,respectively.Conversely,E2F4 might repress TTK expression at the transcriptional level.

Fifty-seven overlapping DEmiRNAs (36 up-and 21 down-regulated)were identified and obtained from datasets GSE114110 and GSE43732(Table S3).A miRNAs-KGs regulatory network comprising 7 KGs,16 overlapping DEmiRNAs (13 up-and 3 down-regulated),and 48 interactions was subsequently constructed (Figure 6C).Survival analysis indicated that the high expression of hsa-miR-483 was negatively correlated with the overall survival of patients with ESCC(P＜0.0001,Figure 6D).Other miRNAs (hsa-miR-183,hsa-miR-145,hsa-miR-139,hsa-miR-708,hsa-miR-7,hsa-miR-96,etc.) presented similar prognostic outcomes in patients with ESCC,although there was no statistical difference (Figure S8).

Targeted drugs analysis

A total of 23 potential anti-ESCC compounds targeting KGs were identified,and most small molecular drugs (20/23) might target TOP2A,AURKA,AURKB,and TTK for the treatment of ESCC (Table 1).Co-inhibitors of AURKA and AURKB included AMG-900,AT-9283,barasertib,danusertib,GSK-1070916,MK-5108,PF-03814735,and SNS-314.

Table 1 Candidate drugs targeting key genes (KGs) in esophageal squamous cell carcinoma (ESCC)

Discussion

Although numerous studies on the molecular mechanism,early diagnosis,and clinical treatment of ESCC have been performed in the past few years,the morbidity and mortality of ESCC are still increasing globally [3].The main reason for the above may be the false positive or negative rates of a single independent analysis.In this study,we used four independent microarray datasets from different reports in the GEO database to obtain more accurate and reliable information.Four hundred and sixty-two overlapping DEGs (212 upand 250 down-regulated) were identified for in-depth exploration of the biological mechanisms and KGs of ESCC (Figure 1E and F).Cell cycle,transcriptional mis-regulation in cancer,ECM-receptor interaction,DNA replication,protein digestion and absorption,proteoglycans in cancer,and biosynthesis of unsaturated fatty acids have been previously shown to be involved in ESCA [27],and we found that they may also participate in the carcinogenesis and progression of ESCC(Figure 5B).Numerous publications have pointed out that abnormal DNA replication and cell cycle are important signals for the initiation and progression of various cancers,including ESCC[28].ECM is an essential ingredient of the cancer cell niche and the first barrier against tumor invasion and metastasis [29].Consistent with previous studies,overlapping DEGs associated with the collagen family in this work (e.g.,COL1A1,COL1A2,COL4A2,and COL4A1)were enriched in the ECM-receptor interaction pathway (Figure 5B),which undoubtedly exacerbated malignancy growth,angiogenesis,and invasion [30].Remodeling of energy metabolism is a symptom of malignancy.Therefore,it is inevitable that ESCC cells create a more suitable environment for growth and invasion by activating pathways related to substance synthesis and energy metabolism,such as protein digestion and absorption,proteoglycans in cancer,and biosynthesis of unsaturated fatty acids [31].

Figure 6 Genetic alterations analysis,and exploration of Transcription factors (TFs) and miRNAs regulatory networks of nine key genes(KGs).(A) Genetic alterations evaluation of KGs.(B) TFs-KGs regulatory network.Red rhomboids: KGs.Red triangles: up-regulated TFs in ESCC tissues compared with normal tissues.Gray triangles: TFs with no significant expression changes in ESCC tissues compared with normal tissues.Delta-shaped arrows:activation of KGs by TFs.Half circular arrows:repression of KGs by TFs.(C)miRNAs-KGs regulatory network.Red rhomboids:KGs.Red rectangles: up-regulated DEmiRNAs.Blue rectangles: down-regulated DEmiRNAs.(D) Prognostic analysis of hsa-miR-483.TFs,transcription factors;KGs,key genes.

NDC80,BUB1,TOP2A,AURKA,AURKB,TTK,UBE2C,TPX2,and BUB1B were identified as KGs in ESCC using multiple machine algorithms,and their expression changes in diseased tissues were validated using a large number of ESCC-related data recorded in the Oncomine and TCGA databases (Figure 2 and Figure S2).The expression of TOP2A and AURKA in ESCC tissues was significantly higher compared to normal tissues at the translational level (Figure 3A-D).The methylation levels of BUB1,UBE2C,and BUB1B were higher in ESCC than those in normal samples (Figure S3).We also evaluated the inhibitory effect of KGs knockout on the survival and growth of 19 ESCC cell lines using genome-wide CRISPR-Cas9 dropout screens (Figure 4).In addition,nine KGs were significantly associated with immune cell infiltration,and they were also strongly positively correlated with each other (Figure S5 and S6),demonstrating that these KGs might be functional complexes that play a vital coordinating role in manipulating the immune microenvironment of ESCC.

We performed a literature review of nine promising biomarkers for ESCC.NDC80 is a mitotic protein that interacts with other proteins to regulate the cell cycle.Studies have shown that NDC80 is over-expressed in colorectal cancer and clear cell renal cell carcinoma[32,33].Based on current findings,we speculate that NDC80 plays a potentially critical role in the carcinogenesis of ESCC.BUB1 and BUB1B have similar functions in the cell cycle.Currently,the role of BUB1 as an oncogene in human cancers including gastric cancer and liver cancer has been observed [34].The clinical value of BUB1B as a potential target for the diagnosis,treatment,and prognosis of hepatocellular carcinoma was previously demonstrated [35].TOP2A is considered a potential biomarker for the diagnosis and prognosis of patients with malignant tumors [36].Both AURKA and AURKB are mitotic regulators,which participate in the regulation of chromosome arrangement and separation during mitosis and meiosis.The over-expression of AURKA is related to the differentiation degree,invasive ability,distant lymph node metastasis,and poor prognosis of ESCC [7].AURKB has been reported to affect the occurrence of ESCC[37].The above studies have pointed out the application direction of AURKA and AURKB,although there is no clear indication that they can be used as clinical diagnosis or treatment targets for ESCC.TTK is a core unit of the spindle assembly checkpoint.Several reports indicate that the abnormally expressed TTK gene is a potential prognostic and/or therapeutic target for human malignancies such as thyroid cancer,breast cancer,liver cancer,and glioblastoma [38].Data from the HPA database showed that UBE2C protein was over-expressed in normal esophageal squamous epithelial cells and moderately expressed in ESCC tissues.Conversely,Palumbo et al.found that UBE2C protein was highly expressed in ESCC samples,and its silencing reduced the malignant phenotype of ESCC cells [39].In recent years,increasing reports have found that high expression of TPX2 further leads to more active proliferation of tumor cells [40].Therefore,it is believed that KGs may be potential oncogenes or driver genes of ESCC,and suppressing their expression can directly or indirectly inhibit the growth and invasion of ESCC cells.

During the exploration of the TFs-KGs network,we identified 12 TFs targeting KGs,of which E2F1,E2F3,E2F4,and MED1 had the most target KGs (Figure 6B).Numerous reports indicate that E2F family members are essential for coordinating cell cycle progression [41].E2F1 is a transcriptional activator of AURKA and AURKB,and its over-expression is a vital indicator of tumor progression and prognosis in patients with ESCC [42].Furthermore,AURKA over-expression increases the stability and transcriptional activity of E2F1 protein in cancer samples by inhibiting proteasome-dependent protein degradation[43].The current study showed that E2F1,as an activator of AURKA,AURKB,and TOP2A genes,aggravated the deterioration of ESCC.E2F4 had a potential reversal effect on ESCC by inhibiting the expression of TTK.Elevated MED1 is a key molecular event related to the occurrence of prostate cancer[44].Similarly,we found that MED1 was abnormally up-regulated in ESCC tissues and might have an activating function for AURKA.

Candidate miRNAs,as oncogenes or tumor suppressors,play important roles in the detection,treatment,and prognosis of malignant tumors [45].Jiang et al.validated that miR-29a-5p could down-regulate the abnormal TPX2 gene,thereby inhibiting the invasion and proliferation of endometrial cancer-derived cells and promoting apoptosis [46].We explored the potential regulatory relationships between KGs and miRNAs to identify miRNA biomarkers related to ESCC (Figure 6C).Consistent with previous studies,we found that hsa-miR-183-5p,hsa-miR-21-5p,hsa-miR-483-5p and hsa-miR-196a-5p had potential diagnostic and/or prognostic biomarker roles for ESCC [47-49].As far as we know,other dysregulated DEmiRNAs might also directly or indirectly affect the etiology and progression of ESCC,and more experiments were needed to support these hypotheses.

Studies have proved that alisertib can inhibit gastrointestinal cancer and glioblastoma by inhibiting the expression of AURKA [50,51].Danusertib is also an effective pan-Aurora kinase inhibitor,which results in cell cycle arrest and programmed cell death in tumor cells[52].Dexrazoxane,an inhibitor of TOP2A,inhibits cancer progression by inducing DNA fragmentation and DNA damage in cancer cells[53].As far as we know,the above drugs are rarely used in the treatment of ESCC,so the data may provide new clues for targeted therapy of ESCC.We speculate that the remaining five KGs may be novel therapeutic targets for ESCC in the future,although no drug information related to them was obtained in this study.

Conclusions

In conclusion,we identified nine KGs(NDC80,BUB1,TOP2A,AURKA,AURKB,TTK,UBE2C,TPX2,and BUB1B)and some significant TFs and miRNAs from large-scale ESCC-related expression data,which may be promising biomarkers for the diagnosis,treatment,and prognosis of ESCC.We also screened 23 potential anti-neoplastic drugs.These data provide more useful diagnostic and therapeutic strategies for clinical research in ESCC.