Soft Tissue Ewing Sarcoma Cell Drug Resistance Revisited: A Systems Biology Approach

Ewing sarcoma is a rare type of cancer that develops in the bones and soft tissues. Drug therapy represents an extensively used modality for the treatment of sarcomas. However, cancer cells tend to develop resistance to antineoplastic agents, thereby posing a major barrier in treatment effectiveness. Thus, there is a need to uncover the molecular mechanisms underlying chemoresistance in sarcomas and, hence, to enhance the anticancer treatment outcome. In this study, a differential gene expression analysis was conducted on high-throughput transcriptomic data of chemoresistant versus chemoresponsive Ewing sarcoma cells. By applying functional enrichment analysis and protein–protein interactions on the differentially expressed genes and their corresponding products, we uncovered genes with a hub role in drug resistance. Granted that non-coding RNA epigenetic regulators play a pivotal role in chemotherapy by targeting genes associated with drug response, we investigated the non-coding RNA molecules that potentially regulate the expression of the detected chemoresistance genes. Of particular importance, some chemoresistance-relevant genes were associated with the autonomic nervous system, suggesting the involvement of the latter in the drug response. The findings of this study could be taken into consideration in the clinical setting for the accurate assessment of drug response in sarcoma patients and the application of tailored therapeutic strategies.


Introduction
Sarcomas represent one-fifth of pediatric cancers and 1% of adult solid malignant tumors. They rise most frequently from connective tissues and only in one-tenth of cases from osseus tissues. Based on their tissue of origin, there are two main classes of sarcomas. The pathogenesis of sarcomas is still unclear, as amid established risk factors reside genetic and epigenetic causes. Literature evidence includes smoking, age, gestational age and weight, parental and maternal health status, occupational exposure to drugs, chemicals, etc. Unfortunately, the survival prognosis of sarcomas is rather poor [1].
At least seventy types of sarcomas have been described so far [2]. Ewing sarcoma is a rare type that develops in bones (skull, spine, pelvis, chest, legs, arms, feet) and soft tissues (head and neck, legs, retroperitoneum) equally frequently and within 13-20 years [3]. More specifically, the Ewing sarcoma family of tumors includes neoplasms localized in the chest wall; this type of peripheral primitive neuroectodermal cancer has been termed an Askin

High-Throughput Gene Expression Data
The publicly available repository NCBI GEO (Gene Expression Omnibus) DataSets (https://www.ncbi.nlm.nih.gov/gds/; accessed on 23 October 2022) [32,33] was thoroughly searched for gene expression datasets related to sarcomas and drug treatment using the keywords: "drug" AND "sarcoma" AND ("human" or "homo sapiens"). The criteria for selecting datasets were: (i) gene expression data from treated and untreated human sarcoma tissues or cell lines, (ii) more than 5000 genes were included in the transcriptomic dataset.
In this way, one eligible RNA-Seq dataset was obtained. The GEO series GSE118871 [34] (Table S1) includes the genome-wide gene expression of Ewing sarcoma A673 cells treated with SP-2509. This dataset contains cell lines both responsive and resistant to SP-2509. The drug-resistant cell lines (herein referred to as "resistant") were established by exposing the corresponding parental A673 cells (untreated) to increasing concentrations of SP-2509 over a 7 month period (in increments of 100 nM) or 48 h (2 uM). Both the resistant cells and those responsive to drug treatment (hereinafter called "responsive") were assessed by Pishas et al. [34] based on experimental cell viability assays; the resistant cells demonstrated a significantly increased viability as compared to parental A673 cells treated with SP-2509, whereas the responsive cells displayed a reduced viability following drug treatment (Table S1). The Illumina HiSeq 2000 (Homo sapiens) GPL11154 platform was employed.
To identify DEGs between resistant versus responsive, and responsive versus parental Ewing sarcoma cells, the EdgeR package version 3.32.0 [39] of the R statistical computation environment v.3.6.1 (https://www.r-project.org; accessed on 28 October 2022) was used. The negative binomial (NB) distribution was used to model the RNA-Seq read counts per gene per sample in EdgeR. Then, the estimating dispersion was calculated with the estimateDisp function. Differential expression analysis between the two RNA-Seq groups was performed using the exactTest function of the EdgeR package v3.32.0. For determining statistically significant DEGs, the cutoff for the absolute log 2 -fold change (FC) was set at two (|log 2 FC≥2|), and the Benjamini and Hochberg (BH)-adjusted p-value [40] was ≤0.05.
The official HUGO Gene Nomenclature Committee (HGNC) [41] gene symbols and gene names were used.

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA), a method to identify biological terms that are enriched in a large gene set, was conducted to functionally annotate the drug resistanceassociated genes detected in this study. To this end, the 'resistant' DEGs were provided as input to WebGestalt (WEB-based GEne SeT AnaLysis Toolkit) [42,43] to identify statistically significant over-represented Gene Ontology (GO) terms; the non-redundant Biological Process subontology of GO was selected, the threshold for the BH-corrected p-value [40] was set at 10 −3 ; the hypergeometric distribution was applied.

Functional Association Network
The associations among the 'resistant' genes/proteins under study were investigated and visualized with the usage of the STRING (Search Tool for Retrieval of Interacting Genes/Proteins) v11.5 [44], a database of both experimentally derived or predicted, direct or indirect, association data among genes/proteins extracted from diverse resources. A relatively high confidence score (≥0.6) for displaying interactions was chosen. The associations were further investigated, analyzed, and visualized through the open-source platform Cytoscape (v.3.8.2) (https://cytoscape.org/; accessed on 11 December 2022) [45].

Epigenetic Regulators of Chemoresistant Genes
The ncRNAs-i.e., miRNAs and lncRNAs-likely regulating the 'resistant' DEGs under study were investigated by employing state-of-the-art software tools.
The DIANA-LncBase v.3 (https://diana.e-ce.uth.gr/lncbasev3; accessed on 24 January 2022) [51], a comprehensive database of experimentally supported miRNA targets on noncoding transcripts, was employed to identify those lncRNAs that potentially interact with the miRNAs detected in the previous step. To increase the accuracy of our analysis, only the miRNA-lncRNA interactions with high confidence detected in cancer/malignant cell types were considered.
An in-house Python script was used to retrieve miRNA-mRNA and miRNA-lncRNA interactions from the respective databases.

Results
The overall procedure followed in this study is outlined in Figure 1.

Epigenetic Regulators of Chemoresistant Genes
The ncRNAs-i.e., miRNAs and lncRNAs-likely regulating the 'resistant' DEGs under study were investigated by employing state-of-the-art software tools.
The DIANA-LncBase v.3 (h ps://diana.e-ce.uth.gr/lncbasev3; accessed on 24 January 2022) [51], a comprehensive database of experimentally supported miRNA targets on noncoding transcripts, was employed to identify those lncRNAs that potentially interact with the miRNAs detected in the previous step. To increase the accuracy of our analysis, only the miRNA-lncRNA interactions with high con dence detected in cancer/malignant cell types were considered.
An in-house Python script was used to retrieve miRNA-mRNA and miRNA-lncRNA interactions from the respective databases.

Results
The overall procedure followed in this study is outlined in Figure 1.  Graphical illustration of the overall methodology of this study. Genes differentially expressed between the chemoresistant and chemoresponsive, as well as the chemoresponsive versus parental Ewing sarcoma cells were detected. Functional annotation analysis of the genes upregulated in the chemoresistant cells and their protein products was performed. Subsequently, a ceRNA was generated from the drug resistance genes, their corresponding regulating miRNAs, and the lncRNAs that act as sponges of these miRNAs.

Gene Expression Profiles of Ewing Sarcoma Drug Resistance
The differentially expressed genes (DEGs) detected between the chemoresistant and chemoresponsive, as well as the chemoresponsive versus the untreated parental Ewing sarcoma cells [34], were 1601 and 1196, respectively (Table S1). Among the genes found upregulated in the drug resistant cells were those encoding "drug resistance-associated membrane proteins" or "DRAMPs", which affected the transport of drugs across membranes. The ABCA2/8/9 genes belong to the broad ATP-binding cassette (ABC) transporter superfamily, the members of which pump drug molecules out of the cell, thereby decreasing the net accumulation of drugs within cancer cells. ABCA2 and ABCA9 were down-regulated in the responsive cells. The expression of a great number of genes coding for another class of DRAMPs, the solute carrier (SLC) transporters, which interfere with the translocation of drug molecules across membranes relying on physicochemical processes [52,53], was also dysregulated (Table S1). A key player in angiogenic signaling in cancer, the vascular epithelial growth factor, VEGFB, was overexpressed in resistant cells (Table S1). Treatment with specific VEGF inhibitors results in transient tumor vascular normalization and increased cancer cell response to chemotherapeutic drugs [54]. An increased expression of ALDH1A3 and ALDH1B1-encoding aldehyde dehydrogenases (ALDH), detoxifying enzymes that catalyze the oxidation of intracellular aldehydes-was detected in the chemoresistant cell lines; ALDHs have been proposed as biomarkers of chemoresistance [55]. ALDH1B1 was also under-expressed in the chemoresponsive cells.
Of note, non-common DEGs (with the same direction, either up-or down-regulated) were detected when the resistance versus responsiveness and responsiveness versus parental (untreated) cells were compared (Table S1). These findings indicate that our results were biologically meaningful, as different molecular factors and mechanisms were implicated in these two phenomena.

Functional Annotation Analysis
A total of 641 'chemoresistant' genes products were found to be implicated in known biological processes ( Figure 2 and Table S2). The protein products of 369 (out of 641) genes formed a highly interconnected network ( Figure 3 and Table S1), suggesting a functional association among these genes leading to their drug resistance effect. One of the overrepresented biological processes was 'drug transport', which included 27 genes, 22 of which were up-regulated in the drug resistant cells ( Figure 2 and Table S1). Given that an increased drug transport activity plays a critical role in drug resistance [7,56] and the bioentities that participate in networks are of a higher biological significance [57], we selected these eight over-expressed genes (DRD1, DRD2, GHRL, KCNA2, SLC7A10, SLC25A13, STRA6, SYT2), the protein products of which participate in the constructed network ( Figure 3). Of note, the genes DRD2, KCNA2, SLC7A10, STRA6, and SYT2 were under-expressed in the chemoresponsive versus untreated parental cells (Table S1).

Epigenetic Regulatory Network of Drug Resistance-Relevant Genes
To gain a glimpse into the post-transcriptional regulatory mechanism(s) of the eight SP-2509-resistant genes, a ceRNA network was constructed, consisting of the miRNAs that regulate the expression of the drug resistance genes, and the lncRNAs, capable of acting as molecular sponges of these miRNAs. For instance, several studies suggest that the downregulation of cancer-relevant miRNAs through sponging alleviates the suppression of downstream mRNAs, affecting different aspects of carcinogenesis [58][59][60].
The miRNAs that potentially target the eight genes in Figure 3 were predicted using different methods. The computational methods for miRNA/mRNA target prediction usually depend on sequence-based features, thermodynamic stability, evolutionary information, or probabilistic models, etc. [61][62][63]. In our study, we used four tools based on different algorithms, so as to extract the maximum possible information. To enhance the prediction accuracy, only those miRNA-target gene interactions predicted by more than three tools were included in the study. Moreover, to obtain more robust results, miRNAs targeting more than two genes were selected for analysis. Collectively, 30 miRNAs were found to target more than 2 genes, and, conversely, 5 genes were targeted by more than 2 miRNAs (DRD1/2, SLC25A13, STRA6, SYT2), suggesting possible co-regulation at the post-transcriptional level. Of these miRNAs, based on our stringent selection criteria, 8 (hsa-miR-149-5p, hsa-miR-29a/b/c-3p, hsa-miR-330-5p, hsa-miR-501-3p, hsa-miR-760 and hsa-miR-766-3p) interacted with the lncRNAs in DIANA-LncBase (Table S2). We retained only the lncRNAs that were upregulated in the drug resistant cells-that is, co-upregulated with the target genes (Tables S1 and S2). Subsequently, a competing endogenous RNA network of the selected six lncRNAs (EDRF1-DT, HAGLR, LINC00997, LOXL1-AS1, SRRM2-AS1, TMPO-AS1), their sponged miRNAs, and the genes targeted by the miRNAs was constructed; however, the mirnas that interact with the DRD2 were not sponged by any lncRNAs (Figure 4).

Discussion
In this work, we investigated the molecular determinants of the Ewing sarcoma cell's line A673 response to diverse anti-neoplastic agents through an in silico approach. To this end, a relevant, publicly available gene expression dataset was used, in which Pishas and colleagues [34] treated Ewing sarcoma cells with increasing concentrations of SP-2509 to investigate the mechanisms underlying the resistance to SP-2509, a small molecular reversible inhibitor of LSD1 [76], in sarcoma. The lysine-specific demethylase 1 (LSD1), also known as KDM1A, demethylates histone 3 lysine 4 (H3K4), thereby acting either as a transcriptional co-repressor or as a transcriptional co-activator by catalyzing the demethylation of histone 3 lysine 9 (H3K9) [77,78]. LSD1 has been demonstrated to contribute to carcinogenesis via chromatin modification, as it promotes or represses the transcription of oncogenes or tumor suppressor genes, respectively [79][80][81][82]. The majority (76%) of the genes found to be differentially expressed between the SP-2509 resistant and responsive Ewing sarcoma cells were up-regulated (Table S1), suggesting that the overexpressed The three members of the mir-29 family, miR-29a/b/c, were the ones mostly targeted by the lncRNAs (Figure 4). In particular, miR-29s have been reported to act primarily as tumor suppressors in numerous cancers through the upregulation of tumor suppressor genes or/and downregulation of oncogenes; they can, thus, elicit apoptosis and inhibit invasion and proliferation of cancer cells [64,65]. In a similar manner, they suppressed the activity of the SP-2509 resistance genes, resulting in the Ewing sarcoma cells responding to drugs. Conversely, the lncRNAs identified in this study could act as drivers of chemoresistance in sarcomas through their sponging role regarding miRNAs (Figure 4). In agreement with this, accumulating evidence suggests that LOXL1-AS1 [66][67][68][69], HAGLR [70,71], LINC00997 [72], and TMPO-AS1 [73][74][75] contribute to carcinogenesis by acting as molecular sponges of miRNAs.

Discussion
In this work, we investigated the molecular determinants of the Ewing sarcoma cell's line A673 response to diverse anti-neoplastic agents through an in silico approach. To this end, a relevant, publicly available gene expression dataset was used, in which Pishas and colleagues [34] treated Ewing sarcoma cells with increasing concentrations of SP-2509 to investigate the mechanisms underlying the resistance to SP-2509, a small molecular reversible inhibitor of LSD1 [76], in sarcoma. The lysine-specific demethylase 1 (LSD1), also known as KDM1A, demethylates histone 3 lysine 4 (H3K4), thereby acting either as a transcriptional co-repressor or as a transcriptional co-activator by catalyzing the demethylation of histone 3 lysine 9 (H3K9) [77,78]. LSD1 has been demonstrated to contribute to carcinogenesis via chromatin modification, as it promotes or represses the transcription of oncogenes or tumor suppressor genes, respectively [79][80][81][82]. The majority (76%) of the genes found to be differentially expressed between the SP-2509 resistant and responsive Ewing sarcoma cells were up-regulated (Table S1), suggesting that the overexpressed genes contributed the most to drug resistance. Of those, eight genes were implicated in drug transport and their corresponding protein products were components of a highly connected network (Figures 2 and 3). These genes/proteins included the signaling receptor and transporter of retinol STRA6, the synaptotagmin 2, the potassium voltage-gated channel KCNA2, as well as the receptor ligands ghrelin and obestatin prepropeptide GHRL. The dopamine receptors DRD1 and DRD2, a class of G protein-coupled receptors, activated adenylyl cyclase to convert ATP into cyclic AMP (cAMP), a second messenger, and were implicated in the dopaminergic regulation of essential neurophysiological processes [83,84]. Two solute carrier (SLC) transporters, SLC25A13 and SLC7A10, were also upregulated. SLC transporters facilitated the influx of drug molecules into cells [85]. However, it has also been suggested that several SLC transporters can mediate bi-directional transport (i.e., both influx and efflux) [86].

Autonomic Nervous System and Drug Resistance
The autonomic nervous system (ANS) regulates many bodily functions, participating in a major fashion in the maintenance of homeostasis and the body's response to stress, and may, thus, influence carcinogenesis. On the other hand, members of the dopamine receptor (DR) family are upregulated in diverse cancers. A positive correlation was observed between a reduced cancer risk and neurodegenerative disorders, such as Parkinson's disease or schizophrenia, where DR-targeting drugs are administered. Moreover, DR antagonists have displayed anticancer efficacy [87,88]. Amid the nodes involved in our network, SYT2 has been suggested to play a regulatory role in synaptic vesicle trafficking and in promoting metastasis in ovarian cancer [89]. In addition, polymorphisms in the gene encoding the neuropeptide GHRL have been associated with non-Hodgkin lymphoma [90]. Three more of the identified major hubs in our constructed interactome: (1) the ion transporter KCNA2, normally expressed in the central nervous system (CNS) [91]; (2) the STRA6, the expression of which was reported in the differentiating nervous system [92]; and (3), the SLC7A10, also called the astrocytic transporter (Asc-1), which has been proposed as a primary driver of the D-serine uptake in the CNS [93]. Taken together, the above findings suggested a link between drug resistance and nervous system function.
Of note, in a study by Gaynes et al., the CNS niche was shown to enhance chemoresistance in leukemia [94]. Additionally, denervation enhanced the effectiveness of chemotherapy in gastric cancer [95]. Finally, Logotheti and colleagues had suggested that genes implicated in neuronal function are activated in cancer cells [96]. As little is known on the potential effects of the ANS on cancer chemoresistance, we compared the up-regulated drug resistance-associated molecules identified herein (Table S1) with murine knockout genes associated with abnormalities in neuronal development, as reported by Logotheti et al. [96]. Of the 119 molecules found in common (Table 1), 40 were oncogenes, including 2 Ewing sarcoma specific genes-namely, the NK2 Homeobox 2 (NKX2) and FEV transcription factor (FEV) [97][98][99]. It has been reported that the NKX2 expression as a biomarker has a high sensitivity (100%), but a moderate specificity in cytologic specimens [97]. Yet, the Ewing sarcoma specificity increased when combined with CD99 [100]. The NKX2 is selectively expressed in the brain, thyroid gland, parathyroid glands, lungs, skin, and enteric ganglia, from prenatal development to adulthood, and has key functions at the interface of the brain, the endocrine, and the immune systems. It is highly expressed in mature limbic circuits related to context-dependent goal-directed patterns of behavior, social interaction and reproduction, as well as in fear responses [101]. Table 1. Common genes between the SP-2509 drug resistance associated genes and murine genes affecting nervous system development. Associated disease syndromes for each gene may be crosschecked in the GeneCards/disorders platform. Genes associated to neurological disorders are marked in blue, genes associated to tumorogenesis or malignancies are marked in red, whilst genes involved in neural and/or cancer entities are marked in purple.

Gene Symbol
Gene   The Nkx2-1 is indirectly involved in autonomic regulation as a major central target for angiotensinogen produced in the subfornical organ and directed to the hypothalamic paraventricular nucleus, where it is converted to angiotensin II [101]. Angiotensin II in the paraventricular nucleus contributes to the autonomic output and sympathetic nervous system excitation [102] and the so termed "enteric branch" (third autonomic branch beyond the sympathetic and parasympathetic ones) of the ANS, as well [103][104][105][106]. The peculiar behavior of the latter branch was initially described by Sternini (1997) [107]. Its association to obesity was first observed in animal models, followed by a meta-analysis in human bariatric surgery data [108], and then described in a preliminary molecular network involved in the "third type of obesity" [109]. Yet, concerning cancer, although this protein has been proposed as a cytostatic factor in certain cancer types, this has not been confirmed on a large scale [110].
FEV gene rearrangements have been observed in a fraction of Ewing sarcoma patients (3.5%) being linked to axial extraskeletal locations (mainly soft tissue), "older age at diagnosis and aggressive clinical behavior", but retaining the classic Ewing sarcoma image [99].
Most of the rest are molecules associated with sensory (deafness, eye disorders) and motor impairments, cardiovascular diseases (i.e., valvulopathies, Brugada syndrome, hypertension), metabolic disorders, neurodegenerative disorders (Parkinson's, Alzheimer, myoskeletal disabilities), encepalopathies and various syndromes, as shown in Table 1. In this table, some molecules are of unidentified function, meriting further research.
Furthermore, five genes involved in drug transport (DRD1/2, SLC7A10, STRA6, and SYT2) and the SLC25A13 paralog, SLC25A1, (Figure 3), are related to neuronal development (Table 1, italicized), further supporting the emerging role of the ANS in drug resistance in chemotherapy.

Conclusions
By employing an interdisciplinary in silico methodology, we identified genes involved in drug resistance in sarcoma cells. A non-negligible fraction of these genes was also associated with the nervous system function and/or development. Of note, the study of the chemoresistance modification effect of the autonomic nervous system in cancer is an ongoing research field with interesting and unpredictable results. Furthermore, this study suggested that future investigations should be directed towards deciphering the crosstalk between the ceRNA and epigenetic regulation of interconnected proteincoding genes. LncRNAs may be considered potential targets for drug resistance based on their capability to attract and inhibit downstream miRNAs, influencing their ability to suppress chemoresistance-associated mRNAs. Both the protein-protein interaction and ceRNA networks comprise of an orchestrated multi-molecular mechanism, the dissection of which would enhance our understanding of the molecular determinants of cancer cell chemoresistance. Therefore, by disrupting the lncRNAs-miRNA interplay and by repressing the corresponding mRNAs, as well as by manipulating interactions between proteins encoded by chemoresistance-relevant genes, the drug resistance of cancer cells could be attenuated. This information can be utilized in the clinical setting for the design of combinatorial therapies, targeting both proteins and epigenetic regulatory factors, such as miRNAs and lncRNAs.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijerph20136288/s1, Table S1: Gene expression data analyzed in this study. Samples from the GSE118871 transcriptomic dataset; differentially expressed genes (resistant versus responsive, responsive versus parental); Table S2: Functional analysis of the Ewing chemoresistant genes.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data and analysis methodologies are contained in the manuscript. Any additional data requests can be addressed to the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.