Identification of key genes and pathways in Ewing’s sarcoma patients associated with metastasis and poor prognosis

Introduction

Ewing sarcoma (ES) is the second commonest primary malignant bone neoplasm accounting for approximately 25–34% of malignant bone tumors. It is a devastating, poorly differentiated, and high-grade osteolytic disease threatening human health severely.^1,2 Statistically, ES is a rare, aggressive bone neoplasm, and around 2.9 people per million encountered with ES annually worldwide, and predominantly appeared in childhood and adolescence.^3,4 ES mainly occurs in bone and surrounding soft tissue characterized by a highly aggressive small round blue cell malignant neoplasm.^5,6 Accumulating evidence strongly documented that chromosomal translocation comprising the Ewing sarcoma breakpoint region (EWSR) gene on chromosome 22 and a member of the ETS family of transcription factors implicated in the ES pathogenetic process, which could encode tumor-specific fusion protein EWSR/FLI and was a distinct and well-defined phenotype for ES genetical characterization.^7–9 Moreover, previous researches revealed that ~13% of patients with ES harbor rare inactivating variants or mutations in DNA damage repair genes consisting of the same genes that are enriched in hereditary breast cancer (such as BRCA1).¹⁰ Currently, surgical excision, multidrug chemotherapy, and local radiotherapy are the principal strategies for the treatment of ES.¹¹ Unfortunately, although enormous progress has been achieved along with science and technological advancements, the diagnose and treatment are far from being satisfactory. Generally, ES patients accompanied with metastatic, or relapsed features have a dismal outcome,^12,13 and the 5-year survival rate for those patients was <30%.¹⁴

To date, emerging progress concerning ES was acquired involving in nuclear architecture and chromosomal positioning, and the actively transcribed genes were located in euchromatin, which tended to cluster to the center of the nucleus.^7,15 Occasional spontaneous DNA breakage may induce the fusion of non-homologous chromosomes and ultimately result in chromosomal translocation.¹⁶ Nevertheless, the pathological mechanisms of ES could not be fully elucidated due to it is intricate and elusive. Therefore, to deeply illuminate the underlying mechanism and related pathways would conduce to thoroughly overcome this disease. In recent years, microarray technology was a great boost for uncovering the pathological mechanisms of ES and substantially accelerating the research phase. Additionally, bioinformatics analysis including patient diagnostic, therapeutic, and pathological information was employed to advance oncology research and lay the foundation for improving disease prevention, early detection, and treatment. The rapid development of bioinformatics enables us to comprehensively screen out the key genes by using high-throughput microarrays. Metastatic status closely correlated with the prognosis of ES patients involving in multiple genes and signal pathways. Effectively and reasonably finding out the hub genes/pathways utilizing microarray database would be in favor of exploiting new strategies for ES treatment. Interestingly, lots of compelling consequences referring to metastatic drivers in ES had been spotted. For example, highly expressed ERBB4 would facilitate tumor metastasis, invasion and suppress apoptosis by activating the phosphoinositide 3-kinase/AKT and focal adhesion kinase pathways in metastatic ES cells.¹⁷ Similarly, differential transcriptomic analyses between primary and metastases revealed that inactive tyrosine-protein kinase transmembrane receptor ROR1 and the putative Wnt family member 5 A ligand jointly participated in the cell migration in ES.¹⁸

In the study, ES metastasis-related abnormally expressed genes were screened out based on public gene expression omnibus (GEO) datasets. Meanwhile, the pivotal biological pathways regarding ES metastasis were clearly defined. Thereafter, biomarkers associated with ES metastasis including mRNA and miRNA were elucidated by conducting a series of bioinformatics and survival analyses. Accordingly, these results may shed a new light on the underlying mechanisms of ES metastasis derived from the bioinformatics.

Materials and methods

Acquirement of microarray data

Microarray data were obtained from the dataset of GSE17618 and GSE12102 gene expression array separately, including primary samples, metastasis samples, and recurrent samples, based on the GPL570 platform ([HGU133 Plus 2] Affymetrix Human Genome U133 Plus 2.0 Array). GSE63157 dataset comprising Ewing’s sarcoma patient’s prognostic information was acquired. GDS4962 dataset containing gene expression profile of ES cell line with inducible EWSR/FLI knockdown were employed to analyze the impact of ESWR/FLI1 knock down over key genes. Moreover, the miRNA expression profile in primary ES and normal human mesenchymal stromal cells (MSCs) was extracted from GEO dataset GSE80201, for upstream miRNA regulator prediction. All datasets are available in the GEO database (https://www.ncbi.nlm.nih.gov/).

Identification of differentially expressed genes (DEGs)

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), an interactive online tool, was employed to compare different sets of samples in the GEO series.¹⁹ To filtrate DEGs, GEO2R was used to filter DEGs between primary and metastasis ES samples based on the GSE17618 and GSE12102 data sets. Thereafter, microarray data of DEGs were downloaded in text format, and the annotation file for probes of the platform was utilized to convert the probe IDs into gene symbols. Thereafter, DEGs were selected for succeeding analysis according to the following criteria: P-value <0.05. Hierarchical clustering analysis was conducted via Morpheus (https://software.broadinstitute.org/morpheus/), and the plot was visualized by RStudio software version 1.1 (R version 3.5.1), with ggplot2 package.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analyses of DEGs

GO is used as a major bioinformatics tool by providing a systematic language to annotate genes, gene products and sequences.²⁰ KEGG is an integrated database resource related to genomes, biological pathways, diseases and chemicals.²¹ The online tool, Database for Annotation Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) provided comprehensive information for a list of genes depending on the GO database and pathway database. The screened DEGs were analyzed using the DAVID tool. P<0.05 was considered statistically significant. Protein–protein interaction (PPI) network analyses were conducted and visualized using GeneMANIA (www.genemania.org/) with the default parameter.²²

Expression and survival analysis of key genes

To investigate the role of key genes in ES’s progression, the expression level including MSH2, MSH6, RPA2, and RFC2 was detected using GSE63157 datasets, respectively. Furthermore, to determine whether these key genes affect clinical outcomes, the series matrix of GSE63157 that contained ES patient’s prognostic information was employed for survival analysis. Briefly, the patients were divided into two groups (high/low group) according to the expression of key genes, with cutoff of average gene expression level. The prognosis of each group was tested by Kaplan–Meier survival estimators, P-value was corrected by Bonferroni method. The relationship between key genes and prognosis value was further analyzed via pan-cancer TCGA database OncoLnc (http://www.oncolnc.org/download/),²³ Cox Regression was conducted. The results were visualized using RStudio software ggplot2 package. Abbreviates for cancer types in TCGA were provided in Table S7. Furthermore, the survival analyses of key genes in TCGA sarcoma (SARC) dataset was conducted and visualized by Oncolnc.

Prediction of upstream regulators

GEO dataset GDS4962 containing 16 samples of A673 ES cell line with or without inducible EWSR/FLI1 knockdown in time course was employed to explore whether EWSR/FLI1 could regulate the expression level of these key genes. Correlation analyses were conducted using R with corrplot package. Additionally, GEO dataset GSE80201 was downloaded for miRNA analyze, DEMs (differentially expressed miRNAs) were acquired by using the parallel method of DEGs mentioned above. Cytoscape app CyTargetLinker consisting of miRTarBase 7.0, TargetScan release 7.2, miRbase version 21, was utilized to investigate the potential correlation between DEMs and DEGs.²⁴

Results

Identification of DEGs

According to differential expression analysis, a total of 298 (Table S1) and 428 DEGs (Table S2) were screened out in metastasis samples based on GSE17618 and GSE12102 dataset when compared to primary samples, respectively. The results showed that 180 genes were upregulated and 118 genes were downregulated in metastasis sample compared to primary samples obtained from GSE17618 (Figure 1A). Moreover, 252 genes were upregulated and 176 genes were downregulated in metastasis sample compared to primary ones in GSE12102 (Figure 1B). Volcano plot of DEGs was selected and showed in Figure1C and D.

Figure 1 Hierarchical cluster analysis and volcano plot of differentially expressed genes (DEGs) in metastasis samples and primary samples. (A and B) represent the DEGs obtained from GSE12102 and GSE17618 dataset, respectively. Each column represents a sample, and each row represents a single gene. Green color represented metastasis samples, and orange color represented primary samples. The color scale shows the relative genes expression level in certain slide: blue indicates low relative expression levels; red indicates high relative expression levels, white indicates zero (no change). (C and D) Volcano plot of DEGs in metastasis samples compared to primary samples. Red indicates the gene expression was upregulated in metastasis samples compared to primary samples (P<0.01); blue indicates the gene expression was downregulated in metastasis samples compared with primary samples (P<0.01); Black indicates the P-value was >0.01.

GO enrichment analysis

The GO database as a dynamic controlled vocabulary is used to describe the role of gene and protein with three categories information consisting of biological process (BP), cellular component (CC), and molecular function (MF). To in-depth investigate the function of the identified DEGs, GO analyses were performed in two datasets (GSE17618 and GSE12102) using DAVID online tool, respectively. In the BP category, the DEGs were mainly involved in neuromuscular junction development, positive regulation of helicase activity, epidermal cell differentiation. CC was significantly enriched in the basolateral plasma membrane, nucleoplasm, SMN-Sm protein complex. Moreover, the GO MFs analysis particularly enriched in deacetylase activity, poly(A) RNA binding, ATP binding. Enriched Go terms with cutoff value of P<0.05 were showed in Table S3 and S4. Top 20 GO terms obtained from GSE12102 and GSE17618 dataset were showed in Figure 2A and B, respectively.

Figure 2 Bioinformatic analysis of DEGs including GO, KEGG and PPI network. (A and B) Bubble plot for GO enrichment of DEGs in GSE12102 and GSE17618 dataset based on GO analysis. X-axis represents the gene ratio and y-axis descript the enrichment components. The area of the circle is proportional to the number of genes assigned to the term and the color accords with the P-value. (C and D) separately showed the KEGG pathways in GSE12102 and GSE17618 dataset. (E) showed the potential intermolecular interactions of the key genes. The different color lines reflect interaction relationship, respectively. (F) Overview of MMR pathway. Pink represents the differently expressed genes (P<0.01); blue represents undifferentiated expressed genes (P>0.01).Abbreviations: DEGs, differentially expressed genes; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; PPI, protein–protein interaction; MMR, mismatch repair.

EGG pathway functional enrichment analyses and PPI network analysis of the DEGs

KEGG analysis result revealed that 10 and 9 terms were involved in ES metastasis as shown in Figure 2C and D. Intriguingly, mismatch repair (MMR), spliceosome, ABC transporters, cysteine and methionine metabolism were simultaneously enriched in both datasets (Figure 2C and D). The most significantly enriched KEGG pathway was MMR pathway (Tables 1, 2). The key genes including MSH2, MSH6, RPA2, and RFC2 that belong to the MMR pathway, were enriched in GSE12102 and GSE17618 dataset. To investigate the PPI relationship, PPI analyses were conducted by using GeneMANIA. The results suggested that PCNA, RFC family of transcription factors, and PRP family proteins have direct or indirect correlation with the 4 key genes (Figure 2E). Moreover, 4 DEGs including MSH2, MSH6, RFC2, RPA2 were annotated in the MMR pathway as shown in Figure 2F.

Table 1 Significant enriched KEGG pathways of DEGs in ES dataset GSE12102

Table 2 Significant enriched KEGG pathways of DEGs in ES dataset GSE17618

The expression and survival analysis of key genes

The expression level and the prognosis value of the key genes on ES patients with/without metastasis were evaluated using GEO datasets and TCGA database. The results showed that the mRNA expression level of MSH2, MSH6, RPA2, and RFC2 was increased in metastasis samples when compared with that of in primary ones (P<0.05, Figure 3A–H). And, the MSH2 expression level was slightly downregulated in recurrence samples compared to metastasis samples (P<0.05, Figure 3A). Meanwhile, GSE63157, and TCGA databases were employed to investigate the prognostic significance of these genes. The result suggested that high expression of MSH2, MSH6, RPA2, and RFC2 was associated with poor overall and event-free survival of ES patients (Figure 4A–H). To investigate the prognosis value of the key genes in pan-cancer TCGA study, Oncolnc.org was applied and COX regression was conducted. The result showed that the expression of the key genes was negatively correlated with the prognosis of SARC, LIHC, LGG, and KIRP patients, and positively correlated with CESC patients (Figure 5A). The full table was presented as  Table S6. Thereafter, survival analyses were conducted in TCGA SARC datasets, the results suggested that MSH2, MSH6, and RPA2 have a significant negative-correlation with SARC patient overall survival (Figure 5B, C, E, P<0.05).

Figure 3 The key genes expression in dataset GSE12102 and GSE 17618. (A–D) showed the expression of key genes including MSH2, MSH6, RPA2, and RFC2 obtained from GSE12102 dataset. (E–H) showed the expression of key genes obtained from GSE17618 dataset. Primary represents primary samples; metastasis represents metastasis samples; recurrence represents recurrence samples. The blue dot represents the gene expression of primary samples; the red square represents the gene expression of metastasis samples; green triangle represents the gene expression of recurrence samples, respectively. *P<0.05; ** P<0.01.

Figure 4 Kaplan–Meier analysis of the key genes in patients with Ewing sarcoma (ES). (A–D) Relationship between MSH2, MSH6, RPA2, and RFC2 expression and the survival of ES patients analyzed using the GSE63157 dataset. (E–H) Relationship between MSH2, MSH6, RPA2, and RFC2 expression and the survival of ES patients analyzed with the GSE17618 dataset. Red curve represents the high expression of the key genes, and the blue curve represents the low expression of the key genes.

Figure 5 Cox-regression and survival analyses of key genes in TCGA database. (A) Cox-regression result of the key genes in TCGA pan-cancer datasets. X axis represents cox coefficient value, positive value indicates high gene expression predicts poor prognoses of patients, negative value represents high gene expression predicts good prognoses; y axis represents cancer types in TCGA database. Different shapes represent different genes, shape size represents the mean expression level of target genes in specific cancer. The color accords with the –log10 (P-value). (B–E) Relationship between MSH2, MSH6, RPA2, and RFC2 expression and the survival of sarcoma patients analyzed using the TCGA SARC dataset. The high expression of the key genes showed in red curve, and the low expression of the key genes showed in blue curve. Abbreviations: TCGA The Cancer Genome Atlas; UCEC, uterine corpus endometrial carcinoma; STAD, stomach adenocarcinoma ; SKCM, cutaneous melanoma; SARC, sarcoma; PRAD, prostate adenocarcinoma; PAAD, pancreatic ductal adenocarcinoma; OV, ovarian serous cystadenocarcinoma; LUSC, lung squamous cell carcinoma; LUAD, lung adenocarcinoma; LIHC, liver hepatocellular carcinoma; LGG, lower Grade Glioma; KIRP, papillary kidney carcinoma ; KIRC, clear cell kidney carcinoma; HNSC, head and neck squamous cell carcinoma; GBM, glioblastoma multiforme; ESCA, esophageal cancer; COAD, colon and rectal adenocarcinoma; CECA, cervical cancer; BRCA, breast cancer; BLCA, urothelial bladder cancer; AML, acute myeloid leukemia.

Evaluation of key genes expression in EWSR/FLI1 knockdown cell and prediction of upstream miRNAs

EWSR/FLI1 expression is closely associated with ES cell viability, and EWSR/FLI1 knockdown promotes the ES cell survival. To further assess the role of key genes including MSH2, MSH6, RPA2, and RFC2 in the prognosis of ES, EWSR/FLI1 knockdown cell was utilized. The result revealed the expression of EWSR and FLI1 was obviously downregulated with time (Figure 6A and B); moreover, the level of MSH2, MSH6, RPA2, and RFC2 decreased first and then slightly increased, and all expression level was significantly lower in EWSR/FLI1 knockdown cell group compared to control group (Figure 6C–F). The correlation was analyzed among different genes, and the results suggested that the expression of EWSR and FLI1 both had a strong association with the four key genes (Figure 6G). Moreover, the ES related miRNA regulation network was predicted using GEO dataset. The results showed that 128 miRNAs were identified as DEMs (Table S5). Network analyses revealed that 24 of them may regulate MSH2, MSH4, RFC2, and RPA2 expression (Figure 6H).

Figure 6 Key genes expression in EWSR/FLI1 knockdown cell and prediction of miRNA Targets. (A–F) Expression of EWSR, FLI1, MSH2, MSH6, RPA2, and RFC2 in EWSR/FLI1 knockdown cell and control group. Black spots represent the control group, and the red spots represent the EWSR/FLI1 knockdown group, respectively. Moreover, the expression level of gene is expressed as different shapes at different time-points in EWSR/FLI1 knockdown group; (G) showed the correlation among different genes. Shades and area of red circle represent the relevancy degree. Arabic numerals represent the person coefficient. The color scales show the correlation. Red represents positive correlation; blue represents negative correlation; (H) showed the miRNA regulation network concerning the key genes predicted by GEO dataset.

Discussion

It has been considered that ES is well characterized in genome level,¹⁰ it is also documented that ES has a very low mutational burden (0.15 mutations/Mb).²⁵ Tumorigenesis is a complex process, multiple levels of gene expression regulation are involved in the process. Transcription level regulation plays a pivotal role in gene expression. However, the potential correlation between transcription level regulation and ES’s metastasis are not clear. Hence, we aimed to further clarify the underlying mechanism regarding gene expression and ES by a bioinformatics approach. In the present study, data obtained from GSE17618 and GSE12102 datasets were analyzed, then GO and KEGG pathway analysis was conducted. The result revealed that the MMR pathway was involved in ES metastasis in both datasets. And, the key genes including MSH2, MSH6, RFC2, and RFC4 were upregulated in metastasis sample compared to the primary sample in both datasets. The key genes are also closely related with the poorer prognosis of ES patients in dataset GSE63157. Meanwhile, TCGA SARC analyses suggested that the key genes upregulation reduced the overall survival in SARC patients. All those results indicated that the key genes play an important role in ES, which may act as prognostic markers for ES.

MMR pathway is highly activated during G1/S-phase,^26,27 in which DNA is duplicated and new histones synthesized. Altered MMR expression indicates a change in proportion of cell in G1/S-phase. Meanwhile, it was reported that overexpression of Msh2 and Msh6 altered the cell cycle distribution in yeast, and result in large percentage of cells in G1, implicating overexpression of MMR genes may promote cancer progression.²⁸ On the other hand, microsatellite instability (MSI) is closely related to DNA MMR deficiency (dMMR). To date, MSI has been identified to affect cancer response to chemotherapies and could regulate PD-1 blockade response in the solid tumor.^29,30 Gryfe et al reported that dMMR tumors had significantly more favorable prognosis,³¹ and this result was supported by other researchers.³² Moreover, Alldinger et al found that dMMR was not prevalent in ES, and concluded that the poor prognosis of ES patients was not associated with the loss of MSH2, MSH6 or MLH1 expression.³³ Similarly, our results further pointed out that the upregulation of MSH2, MSH6 would reduce overall survival in ES patients. Furthermore, the gene rank of MSH2 is 25, whereas MSH6 ranks 704 in SARC TCGA database, which reflect the high impact of MSH2, MSH6 expression on survival of sarcoma patients.

RPA2 is a subunit of replication protein A (RPA) complex, and its dysregulation was observed in various cancers. Studies have revealed the overexpression of RPA2 is associated with poorer survival outcomes and a higher probability of tumor metastasis.^34–36 Also, the RPA2 inhibitor can synergize with DNA-damaging agents to affect chemotherapy response in lung and ovarian cancer cell lines.^37,38 In the present research, RPA2 was upregulated in ES metastasis samples and had a negative correlation with overall survival of ES patients, which indicated that RPA2 may act as a prognosis factor in ES.

Lots of genes participated in this progress of ES metastasis. PPI analyses revealed that multiple key genes directly or indirectly interacted with PCNA. PCNA, as a cofactor of DNA polymerases, encircles DNA and serves as a core element for the recruitment of numerous replication proteins of replication fork.³⁹ PCNA is crucial for Okazaki fragment initiation, and it is loaded onto primer-template junctions by replication factor C (RFC) during DNA replication, after completion of each Okazaki fragment, PCNA is unloaded from DNA, and recycled to promote fidelity of synthesis of subsequent Okazaki fragments.⁴⁰ RFC2 is the only RFC subunit that can independently unload PCNA from DNA. And the expression of RFC2 was significantly elevated in nasopharyngeal cancer,⁴¹ choriocarcinoma,⁴² glioblastoma as well as osteosarcoma.^43,44 RFC2 expression level negatively correlated with ES patient’s survival, which was significantly upregulated in metastasis sample compared to the primary sample. These results suggested that RFC2 may play a key role in the progress of ES metastasis.

EWSR/FLI1 protein regulates multiple target gene expression, which acts as a critical in ES. Heisey et al found that EWSR/FLI1 fusion would increase BCL-2 expression, and thereby induce drug resistant to PARP inhibitors, and BCL-2 and BCL-XL inhibition could significantly reverse the drug resistance in ES.³⁰ FOXO1 is a cancer suppressor in multiple cancer types,^45–47 which significantly inhibited cell proliferation and clone formation ability in ES cell.⁴⁸ However, EWSR/FLI1 fusion could repress FOXO1 expression by binding to the promoter sequence of FOXO1.⁴⁹ Furthermore, knockdown of ESW/FLI1 impairs tumorigenesis of ES in vivo.⁵⁰ In the present study, we identified that four key genes significantly correlated with ES metastasis and overall survival in ES patients. EWSR/FLI knockdown in A673 cell line results in the downregulation of the four key genes, and the expression pattern of EWSR and FLI1 are closely related with the four key genes, which indicated that the upregulation of the four genes may be partially regulated by EWSR/FLI1 fusion. Generally, miRNAs participate in the vast majority of biological regulation. Parafioriti et al reported that miRNAs take part in ES oncogenesis.⁵¹ Our research demonstrated that 23 miRNAs collectively regulated the expression of the four key genes. A pan-cancer Cox regression of the key genes was performed as indicated in Figure 5A, which indicated that the key genes were negatively correlated with the patient’s prognosis. Therefore, the results further documented that the key genes including MSH2, MSH6, RFC2, and RFC4 played a significant role in the progress of ES prognosis.

Our results showed that MMR pathway was significantly correlated with ES metastasis and overall survival. We speculated that the MMR pathway may regulate the proliferation, invasion, and metastasis in ES cells. Whereas lack of biological and clinical experiments to further confirm these results. In the future, relevant experiments comprising clone formation assay, CCK8 assay, wound healing assay, and clinical prospective study will be applied to verify these results.

Conclusion

Our results indicated MMR pathway significantly enricher in ES metastasis samples and may contribute to ES proliferation, invasion and migration. Four key genes of MMR pathway including MSH2 MSH4, RFC2, and RPA2, closely correlated with poor prognosis of ES patients. Further analyses indicated that EWSR/FLI fusion could regulate the key gene expression. Collectively, this research provided a novel insight into the mechanism of ES metastasis and lay a theoretical foundation for prognosis evaluation and chemotherapeutic drug’s selection of ES patients.