Identification of driver genes associated with chemotherapy resistance of Ewing's sarcoma.

BACKGROUND: The aim of this study was to identify the driver genes associated with chemotherapy resistance of Ewing's sarcoma and potential targets for Ewing's sarcoma treatment.
METHODS: Two mRNA microarray datasets, GSE12102 and GSE17679, were downloaded from the Gene Expression Omnibus database, which contain 94 human Ewing's sarcoma samples, including 65 from those who experienced a relapse and 29 from those with no evidence of disease. The differen tially expressed genes (DEGs) were identified using LIMMA package R. Subsequently, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed for DEGs using Database for Annotation, Visualization and Integrated Analysis. The protein-protein interaction network was constructed using Cytoscape software, and module analysis was performed using Molecular Complex Detection.
RESULTS: A total of 206 upregulated DEGs and 141 downregulated DEGs were identified. Upregulated DEGs were primarily enriched in DNA replication, nucleoplasm and protein kinase binding for biological processes, cellular component and molecular functions, respectively. Downregulated DEGs were predominantly involved in receptor clustering, membrane raft, and ligand-dependent nuclear receptor binding. The protein-protein interaction network of DEGs consisted of 150 nodes and 304 interactions. Thirteen hub genes were identified, and biological analysis revealed that these genes were primarily enriched in cell division, cell cycle, and mitosis. Furthermore, based on closeness centrality, betweenness centrality, and degree centrality, the three most significant genes were identified as GAPDH, AURKA, and EHMT2. Furthermore, the significant network module was composed of nine genes. These genes were primarily enriched in mitotic nuclear division, mitotic chromosome condensation, and nucleoplasm.
CONCLUSION: These hub genes, especially GAPDH, AURKA, and EHMT2, may be closely associated with the progression of Ewing's sarcoma chemotherapy resistance, and further experiments are needed for confirmation.

Introduction

Ewing’s sarcoma is a small, round, blue cell tumor and the second most common primary bone cancer in children, adolescents, and young adults.1 The annual incidence of Ewing’s sarcoma is approximately 3 per million in the population.2 The most common primary sites include the long bones (47%), pelvis (26%), chest wall (16%), and spine (6%).2 Ewing’s sarcoma is characterized by a chromosomal translocation involving the EWS gene (EWSR1) on chromosome 22q12 and a member of the ETS transcription factor family. In 85% of patients, the FLI1 gene on chromosome 11q24 is involved, leading to an aberrant EWS–ETS transcription factor.3 Over the last four decades, survival of this disease has improved from 10% with radiotherapy alone to nearly 70% in localized disease with a combination of chemotherapy and local therapy (surgery, radiation therapy, or both).4 Despite the multidisciplinary therapeutic strategies, the treatment of Ewing’s sarcoma has reached a clinical therapeutic plateau. The survival rate of Ewing’s sarcoma has not improved in recent years, and the outcome in some subgroups of patients still remains poor.5 In patients with localized disease, poor histological response or resistance to chemotherapy is associated with poor prognosis.6 Thus, there is an urgent need to identify crucial genes related to chemotherapy resistance of Ewing’s sarcoma and to further develop prognostic biomarkers and novel therapeutic targets.

The current standard-of-care for Ewing’s sarcoma is cytotoxic chemotherapy with vincristine, adriamycin, and cyclophosphamide, alternating with ifosfamide and etoposide.6 However, these agents do not benefit approximately 30% of patients with Ewing’s sarcoma who present with localized disease.7 In addition, those who bear clinically detectable metastases at diagnosis and those who experience local and/or distant recurrent disease still present a low overall survival rate.8,9 Despite multiple efforts to identify new potent therapeutic drugs, there is still a high percentage of tumors that do not respond to treatment. The mechanisms by which Ewing’s sarcoma becomes resistant to chemotherapy are varied and may involve cancer stem cells, proliferative intracellular pathways, and new mutations that allow tumor cells to escape the effect of chemotherapy.10 Therefore, it is imperative to identify prognostic factors to detect chemotherapy-resistant tumors at diagnosis and to develop more individualized treatment regimens.

During the last decades, microarray technology and bioinformatics analysis have been widely used to screen genetic alterations at the genome level, which have helped us identify the differentially expressed genes (DEGs) and functional pathways involved in the chemotherapy resistance of Ewing’s sarcoma. In the present study, two mRNA microarray datasets, GSE12102 and GSE17679, were downloaded from the Gene Expression Omnibus (GEO) database, and a total of 94 human Ewing’s sarcoma samples were selected from these two datasets.11,12 Then, we used LIMMA package R to analyze and obtain DEGs between tissues from patients who experienced a relapse (REL) and from those with no evidence of disease (NED). Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses and protein–protein interaction (PPI) network analyses were performed to elucidate the molecular mechanisms underlying chemotherapy resistance. In conclusion, a total of 347 DEGs and 13 hub genes were identified, and three significant genes were found. These three genes may be crucial in Ewing’s sarcoma resistance to chemotherapy and candidate biomarkers to predict the prognosis of Ewing’s sarcoma.

Methods

Affymetrix microarray data

The mRNA expression profiles of GSE12102 and GSE17679 were downloaded from the public functional genomics data repository GEO (http://www.ncbi.nlm.nih.gov/geo/).11,12 These datasets were obtained on the Affymetrix GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA). A total of 94 human Ewing’s sarcoma samples were included in these two data-sets, and we defined the patients who had tumor progression within 3 years from diagnosis as the REL group and those who did not experience recurrence as the NED group. We then compared the gene expression profiles of 65 tumors from the REL group with 29 tumors from the NED group.

Identification of DEGs

Background correction, quantile normalization, and probe summarization of the raw data in the CEL files were performed using the robust multiarray average algorithm in the affy package of Bioconductor (http://www.bioconductor.org/). The DEGs between REL group samples and NED group samples were identified by the linear models for microarray data (LIMMA, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package. The corresponding P-value of the genes after independent t-tests was defined as the adjusted P-value, and an adjusted P-value <0.01 and |log2FC (fold change)| ≥1 were selected as the cut-off criteria.

KEGG and GO enrichment analyses of DEGs

Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) is an online biological information database that integrates biological data and analysis tools, and it provides comprehensive functional annotation information for genes and proteins for users to extract biological information. GO function and KEGG pathway enrichment analyses of DEGs were performed based on DAVID. P-values <0.05 and counts ≥2 were considered the cut-off criteria.

PPI network construction and module analysis

The PPI network was predicted using Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/). In the present study, the PPI network of DEGs was constructed based on the STRING database and was selected based on a score (median confidence) >0.4. The plug-in Molecular Complex Detection (MCODE) (version 1.4.2) of Cytoscape is an APP for clustering a given network based on topology to find densely connected regions. The network was then visualized by Cytoscape (http://www.cytoscape.org/), and the most significant module in the PPI networks was identified using MCODE. The criteria for selection were as follows: MCODE scores >5, degree cut-off =2, node score cut-off =0.2, max depth =100, and k-score =2. Subsequently, the KEGG and GO analyses for genes in this module were performed using DAVID.

Hub gene selection and analysis

The hub genes were selected with degrees ≥10. The biological process analysis of hub genes was performed and visualized using the Biological Networks Gene Oncology tool (version 3.0.3) plug-in of Cytoscape. To describe the importance of nodes in the PPI network, three methods, including degree centrality, betweenness centrality, and closeness centrality, were used in the present study.

Results

After standardization of the microarray results, a total of 206 upregulated DEGs and 141 downregulated DEGs were obtained from REL samples compared with the NED samples. The hierarchical cluster analysis of the data revealed that the DEGs can be used to accurately distinguish REL samples from NED samples (Figure 1).

Figure 1

Heat map for the DEGs.

Notes: Each row represents a single gene, and each column represents a tissue sample. The gradual color change from green to red represents the change from upregulation to downregulation.

Abbreviation: DEGs, differentially expressed genes.

Enrichment analyses of DEGs

GO functional enrichment analysis showed that upregulated DEGs were primarily enriched in DNA replication, nucleoplasm and protein kinase binding for biological processes, cellular component and molecular functions, respectively (Table 1). The downregulated DEGs were predominantly involved in receptor clustering, membrane raft, and ligand-dependent nuclear receptor binding (Table 2). According to KEGG pathway enrichment analysis, the upregulated DEGs were primarily enriched in pathways, such as antigen processing and presentation, and there was no pathway enrichment in the downregulated DEGs (Table 3).

Table 1

The top five enriched GO terms of the upregulated DEGs

Category	ID	Term	P-value	Gene
BP	GO:0006260	DNA replication	0.001716486	RECQL4, RAD1, TIMELESS, GINS4, KIAA0101, CHTF18, CDC25A, BARD1
	GO:0007067	Mitotic nuclear division	0.001826818	TIMELESS, CCNF, CDCA2, CENPW, AURKA, CENPT, ANLN, CIT, CDC25A
	GO:0051301	Cell division	0.016233622	TIMELESS, NCAPG, CCNF, CDCA2, CENPW, AURKA, CENPT, CDC25A
	GO:0051276	Chromosome organization	0.018439446	CENPW, CENPT, NCAPH2
	GO:0021831	Embryonic olfactory bulb interneuron precursor migration	0.022031365	ARX, RAC1
CC	GO:0005654	Nucleoplasm	2.35E-05	E2F1, HNF1A, AURKAIP1, SYNCRIP, AURKA, SART3, SKAP2, PTMA, CBX5
	GO:0005694	Chromosome	1.54E-04	RECQL4, RAD1, CENPB, CDCA2, NUSAP1, NCAPH2, EHMT2, NSD1
	GO:0005737	Cytoplasm	4.96E-04	HNF1A, TCOF1, LRRC8D, EIF5B, RANGAP1, SKAP2, SART3, PTMA
	GO:0005634	Nucleus	0.002528025	E2F1, HNF1A, AURKAIP1, TCOF1, EIF5B, TRMT10A, AURKA, SART3
	GO:0000776	Kinetochore	0.011939843	CENPW, RANGAP1, CENPT, CBX5
MF	GO:0019901	Protein kinase binding	0.010341955	E2F1, FOXM1, RAC1, TELO2, ELAVL1, AURKA, GPRC5B, CDC25A
	GO:0005515	Protein binding	0.016084872	AURKAIP1, SYNCRIP, RANGAP1, AURKA, RNF216, SART3, CBX5
	GO:0044822	Poly(A) RNA binding	0.016099213	FBRSL1, GTPBP4, FAM98B, STRAP, TCOF1, EIF5B, NUSAP1, ELAVL1
	GO:0003682	Chromatin binding	0.033200105	ARX, ASCL1, CENPB, MTA3, KIAA0101, ATAD2, LRWD1, HDGFRP2, CBX5

Notes: P-value <0.05 is considered statistically significant. There are only four significant GO terms in MF.

Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, gene ontology; MF, molecular function.

Table 2

The enriched GO terms of the downregulated DEGs

Category	ID	Term	P-value	Genes
BP	GO:0043113	Receptor clustering	0.018378	GRIK5, DLG2, CACNA1A
CC	GO:0045121	Membrane raft	0.006816	RET, LTB4R, HPSE, BACE1, ERLIN2, CD226
	GO:0043025	Neuronal cell body	0.009808	CALCA, EFHC1, RET, KCND2, PPARGC1A, CAMK2N1, CACNA1A
	GO:0043235	Receptor complex	0.037073	RET, PLA2R1, INSR, KLRD1
MF	GO:0016922	Ligand-dependent nuclear receptor binding	0.006912	NCOA2, SMARCE1, PPARGC1A
	GO:0030374	Ligand-dependent nuclear receptor transcription	0.037376	SS18, NCOA2, PPARGC1A

Note: P-value <0.05 is considered statistically significant.

Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, gene ontology; MF, molecular function.

Table 3

KEGG pathway enrichment analysis of DEGs

Category	ID	Term	P-value	Gene
KEGG	hsa04612	Antigen processing and presentation	6.02E-04	KIR2DS2, KIR2DS1, KIR2DL1, KIR2DS4, KIR2DL3, KIR2DL2
	hsa04650	Natural killer cell mediated cytotoxicity	8.00E-04	RAC1, KIR2DS2, KIR2DS1, KIR2DL1, KIR2DS4, KIR2DL3, KIR2DL2
	hsa03013	RNA transport	8.96E-04	XPOT, POM121, STRAP, CYFIP2, EIF5B, RANGAP1, POM121C, THOC5

Note: A P-value <0.05 is considered statistically significant.

Abbreviations: DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The PPI network consisted of 150 nodes and 304 interactions (Figure 2). A significant network module was obtained from the PPI network (Figure 3). Strikingly, all the genes in this module were the upregulated DEGs in Ewing’s sarcoma samples. The functional analysis of genes involved in this module was performed using DAVID. The results showed that these DEGs were primarily involved in mitotic nuclear division, mitotic chromosome condensation, and nucleoplasm (Table 4).

Figure 2

PPI network constructed with the up- and downregulated DEGs.

Note: Red nodes represent upregulated genes, and blue nodes represent downregulated DEGs.

Abbreviations: DEGs, differentially expressed genes; PPI, protein–protein interaction.

Figure 3

The significant module identified from the PPI network.

Note: Red nodes represent upregulated genes.

Abbreviation: PPI, protein–protein interaction.

Table 4

Enriched GO terms of DEGs in the module interaction network

Category	ID	Term	P-value	Genes
BP	GO:0007067	Mitotic nuclear division	0.005736	CDCA2, AURKA, ANLN
	GO:0007076	Mitotic chromosome condensation	0.007125	NCAPG, NUSAP1
	GO:0051301	Cell division	0.011163	NCAPG, CDCA2, AURKA
	GO:0045840	Positive regulation of mitotic nuclear division	0.012322	NUSAP1, AURKA
	GO:0000281	Mitotic cytokinesis	0.013736	NUSAP1, ANLN
CC	GO:0005654	Nucleoplasm	2.68E-04	FANCI, FOXM1, KIAA0101, CDCA2, ATAD2, AURKA, ANLN
	GO:0005876	Spindle microtubule	0.019156	NUSAP1, AURKA

Notes: P-value <0.05 was considered statistically significant. There was no MF enrichment in this module.

Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, gene ontology; MF, molecular function.

A total of 13 genes were identified as hub genes with degrees ≥10. The names, abbreviations and functions for these hub genes are shown in Table 5. Biological analysis showed that these hub genes were enriched in cell division, cell cycle and mitosis (Figure 4). Furthermore, betweenness centrality, degree centrality and closeness centrality were used to evaluate the top 15 significant genes in the PPI network (Table 6). We compared the top four genes and found that GAPHD, AUPKA and EHMT2 are the most important genes. These three genes were upregulated DEGs and among the selected hub genes. Thus, we speculated that GAPHD, AUPKA and EHMT2 play a significant role in the chemotherapy resistance of Ewing’s sarcoma.

Table 5

Functional roles of 13 hub genes with degree ≥10

Number	Gene symbol	Full name	Function
1	ATAD2	AAA domain containing 2	May be a transcriptional coactivator of the nuclear receptor ESR1 required to induce the expression of a subset of estradiol target genes, such as CCND1, MYC, and E2F1
2	AURKA	Aurora kinase A	Mitotic serine/threonine kinases that contributes to the regulation of cell cycle progression
3	FOXM1	Forkhead box M1	Transcriptional factor regulating the expression of cell cycle genes essential for DNA replication and mitosis
4	CDC25A	Cell division cycle 25 homologue A	Tyrosine protein phosphatase which functions as a dosage-dependent inducer of mitotic progression
5	RAC1	Ras-related C3 botulinum toxin substrate 1	Rho family, small GTP-binding protein Rac1
6	NUSAP1	Nucleolar and spindle associated protein 1	Microtubule-associated protein with the capacity to bundle and stabilize microtubules
7	NCAPG	Non-SMC condensin I complex, subunit G	Regulatory subunit of the condensin complex, a complex required for conversion of interphase chromatin into mitotic-like condensed chromosomes
8	GAPDH	Glyceraldehyde-3- phosphate dehydrogenase	Has both glyceraldehyde-3-phosphate dehydrogenase and nitrosylase activities, thereby playing a role in glycolysis and nuclear functions, respectively. Participates in nuclear events including transcription, RNA transport, DNA replication, and apoptosis
9	EHMT2	Euchromatic histone-lysine N-methyltransferase 2	Histone methyltransferase that specifically mono- and dimethylates “Lys-9” of histone H3 (H3K9me1 and H3K9me2, respectively) in euchromatin
10	CCNF	Cyclin F	Acting as an inhibitor of centrosome reduplication
11	FANCI	Fanconi anemia, complementation group I	Plays an essential role in the repair of DNA double-strand breaks by homologous recombination and in the repair of interstrand DNA cross-links
12	ANLN	Anillin, actin binding protein	Essential for the structural integrity of the cleavage furrow and for completion of cleavage furrow ingression
13	KIAA0101/PCLAF	PCNA-binding protein	Acts as a regulator of DNA repair during DNA replication

Figure 4

Biological process analysis of the hub genes.

Notes: The biological process analysis of hub genes was constructed using BiNGO. The color depth of the nodes indicates the corrected P-value of ontologies. The size of the nodes indicates the number of genes involved in the ontologies. P<0.01 was considered statistically significant.

Abbreviations: BiNGO, Biological Networks Gene Oncology; GO, Gene Ontology.

Table 6

Top 15 genes evaluated by subgraph centrality, betweenness centrality, and degree centrality in the PPI network

Gene name	Betweenness centrality	Gene name	Degree centrality	Gene name	Closeness centrality
GAPDH	10,420.99	GAPDH	26	GAPDH	115
EHMT2	5,190.157	AURKA	23	EHMT2	40
RAC1	4,263.988	EHMT2	20	CCNF	27
AURKA	3,387.611	NCAPG	17	AURKA	25
CCNF	2,257.71	FOXM1	15	RAC1	20
FOXM1	1,856.824	RAC1	13	NCAPG	13
PPARGC1A	1,481.229	CCNF	13	INSR	11
NCAPG	1,407.738	CDC25A	12	EIF4A2	9
ILF3	1,373.423	NUSAP1	12	SYNCRIP	9
GTPBP4	1,355.005	ATAD2	11	CDC25A	7
SYNCRIP	1,351.062	FANCI	11	HIST1H4D	7
INSR	1,181.45	KIAA0101	11	ILF3	6
HIST1H4D	1,179.832	ANLN	10	TMPO	6
CDC25A	1,087.077	ILF3	9	GTPBP4	6
EIF4A2	1,041.513	GTPBP4	9	CALCA	6

Note: A combined score >0.4 was considered a threshold value for a significant difference.

Abbreviation: PPI, protein–protein interaction.

Discussion

Ewing’s sarcoma is the second most common malignant bone tumor affecting children and adolescents worldwide. Patients with Ewing’s sarcoma are commonly treated with aggressive cytotoxic chemotherapy and local surgery, and the prognosis for patients is poor. In addition, the 5-year survival rate is poor. Although chemotherapy has been used to treat patients since the 1970s,4 35% of the patients undergoing treatment show chemotherapy resistance. Therefore, there is an urgent need to explore the crucial genes and mechanisms of chemotherapy resistance, which would help develop effective diagnostic and treatment strategies.

In the present study, two mRNA microarray datasets were analyzed to obtain DEGs between REL tissues and NED tissues. A total of 273 DEGs were identified among the two datasets, including 206 upregulated DEGs and 141 downregulated DEGs. The upregulated genes were primarily enriched in DNA replication, mitotic nuclear division, and cell division, while the downregulated genes were predominantly enriched in receptor clustering. GO enrichment analysis revealed that changes in the most significant modules were primarily enriched in mitotic nuclear division, mitotic chromosome condensation, and cell division.

Most cytotoxic cancer therapies and radiotherapies can damage DNA or disrupt mitosis to induce cell death in highly proliferating tumor cells. In the present study, the upregulated genes were primarily enriched in DNA replication (RECQL4, RAD1, TIMELESS, GINS4, and KIAA0101), mitotic nuclear division (TIMELESS, CCNF, CDCA2, CENPW, and AURKA), and cell division (TIMELESS, CCNF, CDCA2, CENPW, and AURKA). We hypothesized that these genes may be involved in the progression of chemotherapy resistance in Ewing’s sarcoma by disturbing DNA replication, mitosis, and division of cancer cells. The TIMELESS gene is highly conserved and is involved in cell survival after damage or stress, increased DNA polymerase epsilon activity, maintenance of telomere length, and epithelial cell morphogenesis.13 Changes in the TIMELESS gene or its expression may promote the development of various cancers.14–19 TIMELESS contributes to the progression of breast cancer through activation of MYC.20 TIMELESS forms a complex with PARP1 and plays a role in the DNA damage response.21 One study has reported that overexpression of TIMELESS in NPC cell lines results in resistance to cisplatin-induced apoptosis in vitro and in vivo via the Wnt/β-catenin pathway and downstream gene transcription.22 Another study has reported that TIMELESS knockout sensitizes HCT116 colon cancer cells to doxorubicin.23 In our study, TIMELESS was overexpressed in REL Ewing’s sarcoma tissue. In addition, this gene was enriched in the DNA replication, mitotic nuclear division, and cell division categories. Accordingly, we hypothesized that TIMELESS may be involved in the progress of Ewing’s sarcoma chemotherapy resistance via the DNA replication process.

Furthermore, in the PPI network module, 13 DEGs were selected as hub genes with degrees ≥10. Biological analysis showed that these hub genes were enriched in cell division, cell cycle, and mitosis, indicating that chemotherapy resistance of Ewing’s sarcoma is associated with these processes. We then used betweenness centrality, degree centrality, and closeness centrality to evaluate the most significant genes in the PPI network, and we compared the top four genes and found that GAPHD, AUPKA, and EHMT2 were the most important genes. GAPDH is commonly known as a glycolytic enzyme,24 and is frequently used as a loading control in gene expression and protein studies. However, recent reports have identified some new roles for GAPDH in a diverse range of cellular processes, particularly its roles in the nucleus and in the regulation of cell proliferation.25,26 In the nucleus, GAPDH participates in the transcription of genes involved in antiapoptotic pathways and cell proliferation through its direct interaction with DNA and/or proteins.27–29 GAPDH also plays a role in the regulation of telomere length, and one report has shown that GAPDH protects against telomere shortening induced by the chemotherapeutic agents gemcitabine and doxorubicin.30 Several anticancer therapies, such as structural analogs of nitrogen bases and nucleosides, can disturb the structure and function of nucleic acids, which leads to tumor cell death. These damages can be recognized by different DNA repair systems. One report has shown that GAPDH can interact with proteins that recognize these damages and participate in DNA repair.31 Therefore, we concluded that in Ewing’s sarcoma, GAPDH is involved in DNA repair and plays a role in chemotherapy resistance. In this study, AURKA was enriched in mitotic nuclear division and cell division as well as in the MCODE module. AURKA is a serine/threonine kinase associated with cell cycle maintenance of stemness and regulation of mitotic spindle formation.32,33 Many studies have reported that AURKA overexpression induces centrosome amplification, chromosomal instability, and oncogenic transformation.34–36 In ovarian cancer, esophageal squamous cell carcinoma, breast cancer, and urinary tract urothelial cancer, increased AURKA is correlated with tumourigenesis, clinical aggressiveness, and tumor progression.37–40 Bioinformatics analysis has shown that AURKA is a potential biomarker for non-small-cell lung cancer (NSCLC) and can predict the prognosis of NSCLC patients.41 Upregulation of AURKA induces NSCLC initiation and progression through the LKB1/AMPK signaling pathway.42 In breast cancer, AURKA has been reported to be a possible marker for endocrine resistance.40 Thus, we speculated that AURKA may play an important role in Ewing’s sarcoma chemotherapy resistance and may act as a potential biomarker for prognosis. G9A/EHMT2 is a histone lysine methyltransferase that specifically affects the mono- and dimethylation of Lys9 of histone H3.3–5. Studies have shown that EHMT2 overexpression is correlated with clinicopathological parameters and poor survival in various cancers.43–46 Elevated EHMT2 induces development and progression of various cancers by regulating cancer metabolism, metastasis, cell survival, and response to hypoxia.47–50 In colon cancer-initiating cells, one report has found that EHMT2 mandates self-renewal, chemotherapy resistance, and metastatic abilities by activating the K-RAS/β-catenin pathway.51 Another study has shown that EHMT2 controls colon cancer stem cell phenotype and chemoradioresistance through the PP2A2-RPA axis-mediated DNA damage response.52 In the present study, EHMT2 was an upregulated DEG in Ewing’s sarcoma. Thus, we hypothesized that EHMT2 may play an important role in Ewing’s sarcoma chemotherapy resistance and may act as a potential biomarker of prognosis.

In conclusion, the present study was designed to identify DEGs that may be involved in the chemotherapy resistance of Ewing’s sarcoma. A total of 347 DEGs and 13 hub genes were identified. Among the hub genes, GAPDH, AURKA, and EHMT2 may be highly correlated with chemotherapy resistance and may be diagnostic biomarkers for Ewing’s sarcoma. However, there were certain limitations in this study, and further experiments are needed to validate the expression levels of these genes and their functions.