This article is mentioned in:

Introduction

Ewing sarcoma is a rare soft or bone tissue tumor mainly occurring in the second decade of life and appears to be a slightly heterogeneous group of cancers mainly composed of a fusion of EWSR1 with a variety of ETS transcription factors at the molecular level (1). The specific cells of origin leading to Ewing sarcoma tumors in a true in vivo situation remain undetermined, but there is certain established evidence that Ewing sarcoma arises from mesodermal structures or respectively from neural crest-derived mesenchymal stem cells (MSCs) (2-4). As in numerous cancer entities, studies indicate a different incidence in different human populations (5). For instance, patients with a European ancestry exhibit a 9-fold higher risk compared to African Americans and Asian/Pacific Islanders (0.155 in the European region vs. 0.017 in the Asian/Pacific region and people of African ancestry). In spite of certain progress in the treatment of Ewing sarcoma in recent years, the overall survival rate is between 61 and 74% for patients with a localized disease with deviations in both directions concerning body localization and tumor progression (6). Several therapeutic concepts have been developed up to now (7) and a comprehensive overview of the 'genomic, epigenomic and transcriptomic landscapes' was presented by Sand et al (8).

In contrast to the clinical sphere, where considerable progress has been made in the description, categorization and treatment of Ewing sarcoma, much less basic knowledge has been generated concerning the systems biology of Ewing sarcoma and its progression from the cell of origin to a life-threatening tumor. In 2009, Suvà et al (9) was able to isolate a subpopulation of CD133+ Ewing cells and demonstrate their cancer stem cell (CSC) and MSC properties. Further efforts by Yang et al (10) to isolate Ewing CSC inspired our group to establish and characterize, in an extended way, a similar side population (SP) of Ewing sarcoma cells of the EWSR1-ERG type [cf. fusion gene types (11)]. The SP was comprehensively characterized in several assays, which highlighted their stem- ness character. In conjunction, a differential gene expression (DE) study was performed by our group. The resulting pattern, with a high content of oncogenes, pointed to a plausible alteration in the pathways of AP-1 complex, APC/c-CDC20 and HDAC9, partly explaining the tumor mechanisms in this fully developed sarcoma (12). The biology of the CSCs and their interplay with the main (non-CSC) Ewing cancer population appear to be complex [e.g., (13)], but the molecular patterns highlighted in a previous study by our group indicated rather small but distinct molecular differences between the CSC and non-CSC populations (12).

Describing a DE pattern between an Ewing sarcoma CSC SP and bone marrow-derived MSCs is only one molecular perspective in the complex systems biology of the cell. Another important aspect concerns splice variants, whose contribution to this altered network situation has not been conclusively explored so far, particularly for EWSR1-ERG based Ewing sarcomas.

It is generally accepted that alternative splicing, here in the sense of differential splicing (DS) events have an enormous role in the regulatory response spectrum (14,15) the cellular system has developed, and is essential to cope with a multitude of environmental challenges, and to synchronize multicellular behavior. A deep insight into alternative splicing in cancer is given by Oltean and Bates (16).

Thus, analyzing the splice variants of a whole transcriptome will definitely provide regulatory and response information for the whole transcriptome, which may otherwise remain more limited. For Ewing sarcoma, this information has, to date, been patchy. Patócs et al (17) looked at the breakpoint variants of the Ewing fusion gene itself and observed 9 different variants and even multiple variants in the same tumor sample, also indicating the existence of different response activity levels of one and the same fusion gene. A study by Sand et al (18) focused on hub genes with clinical relevance, such as CXCR4. Downstream target effects of the fusion transcripts on the alternative splicing of further genes such as ARID1A were observed by Selvanathan et al (19). However, up to now, further high-throughput studies adding whole-genome DS information to existing whole-genome DE information are lacking for Ewing sarcoma with the EWSR1-ERG fusion gene.

Studies on alternative splicing may be categorized into those that are mainly interested in screening and mapping the diversity of observable splice variants to form a biological network and those focusing on the splicing process itself (20), although both perspectives intermingle at the network level. The focus of the present study is on the first type.

Consequently, the main objective of the present study was to define the DS information on a whole-transcriptome basis using the CADO-ES1 (CADO) model cell line harboring the fusion gene EWSR1-ERG, a CSC SP thereof, bone marrow-derived MSC primary cell cultures as putative progenitor cells of Ewing sarcoma and further controls. Comparison between MSCs and the SP should enable the detection of similar and dissimilar (stem) cell features. To make the primary splice detection step more robust, rMATS (21) was chosen, a DS tool modeling Bayesian statistics around a differential approach, in contrast to the majority of splice detection tools, which are working in a non-differential way. The performance of rMATS has been convincingly demonstrated in a comparative study from 2017 (22) and also in a newer study from 2019 (23). The downstream analysis of the DS data in the present study focused on, inter alia (i.a.), creating specific DS gene sets, characterizing the expression profiles of wild-type fusion gene partners and annotating the major biological pathways involved in this multi-group design. Finally, these results were compared with established knowledge and an existing DE analysis of the same dataset from 2018 (12) to create an integrated view.

Materials and methods

Samples

Using the Ewing sarcoma cell line CADO-ES1 (CADO), which was established from the malignant pleural effusion of a 19-year-old Japanese female in 1991 (24), the CSC features of an SP were analyzed in comparison to the main Ewing sarcoma population (NSP) and MSCs. The SP was established by fluorescence-activated cell sorting of the cells and extensively characterized by several assays in a previous publication by our group where additional details are provided (12). Furthermore, 8 different human samples with several biological replicates were considered for the present study. A total of 4 samples were based on the CADO cell line (DSMZ-German Collection of Microorganisms and Cell Cultures) carrying the EWSR1-ERG fusion gene: i) SP, three biological replicates; ii) NSP including 1-3% SP cells, three biological replicates; iii) no stain (and no sort) control (nst), three biological replicates; iv) not sorted (but stained) control (nso), three biological replicates; v) primary cell culture of bone marrow-derived MSCs (from two females and two males; age, six to 29 years; bone marrow classified as normal from Ewing patients; March 2017; Department of Pediatric Hematology and Oncology, University Hospital Münster, Münster, Germany), four biological replicates [protocol (25), assessment according to (26)]; vi) two primary Ewing cell cultures (ES): a) EWSR1-ERG (ES1) and b) EWSR1-FLI1 (ES2), with one sample each [DC-ES-6 and -15, respectively (25,27); renamed as MS-EwS-16/-15 (28); two males; age, 12 and 22 years; Ewing sarcoma; March 2017; Department of Pediatric Hematology and Oncology, University Hospital Münster, Münster, Germany]. As a further off-target control, human skin fibroblasts (hFIB) were utilized with four biological replicates [protocol (29), source (30); all four donors are male; age, one to two years; foreskin fibroblasts; March 2017; Centre of Reproductive Medicine and Andrology, University Hospital Münster, Münster, Germany]. Nso means no fluorescence-activated cell sorting of the cells, while nst means that the Hoechst dye to sort the cells was not added to the cells and they were not sorted at all. Both mentioned procedures cause stress for the cells of interest and this may produce false-positive/false-negative DS genes. The respective controls nst and nso should intercept those misleading results. Experiment-related details are mentioned below in the section on experimental design. All patients were included into the multicenter E.U.R.O Ewing 99 (EE99; no. NCT00020566, 12/02/1999) and Ewing 2008 (no. NCT00987636, 01/10/2009) trials.

Sequencing and raw data processing

The dataset was created in-house for a DE study focusing on CSC features of a distinct Ewing sarcoma SP of CADO cells and the data generation was described in detail by Hotfilder et al (12). In brief, for sequencing, 0.5-10×106 cells/ml were used depending on the approach and the applied downstream procedure. The RNA library preparation for the SOLiD5500×l technology was according to a paired-end design. The read length was 75/35 nucleotides, trimmed to 35/35 to retain a high quality for translating color-space into base-space data. The conversion tool is available online in two different flavors (http://complex-systems.uni-muenster.de/tools.html). Approximately 30 million high-quality reads per sample were remaining after preprocessing for the downstream analysis, which denotes an average transcriptome coverage of 35×. The entire dataset E-MTAB-6067 is available online at the EMBL-EBI data repository (https://www.ebi.ac.uk/arrayex-press/experiments/E-MTAB-6067/). An overview of the full data workflow is provided in Fig. 1.

Figure 1 Flow chart providing an overview of what type of data processing was considered. The raw reads were aligned to the human genome (GRCh38) and the resulting BAM files were used for the rMATS analysis. Details of the sample combinations are presented in Fig. 7. Every rMATS result table is based on five basic splice events, either considering splice junction counts only or splice junction and adjacent exon counts. Further processing was performed on the R platform and the R workspace is provided as supplementary information. GO, gene ontology; FDR, false discovery rate; RNAseq, RNA sequencing; rMATS, program for splice variant analysis; DS, differential splicing; PSI, percentage spliced in-DS quality criterion; MXE, mutually exclusive exons; SE, skipped exon; RI, retained intron; A3SS or A5SS, alternative 3′ or 5′ splice site.

EWSR1-ERG fusion gene and wild-type genes

The algorithm of rMATS does not allow for the direct identification of gene fusions; therefore, STAR-Fusion (31) was used to confirm the EWSR1-ERG fusion gene in the CADO samples SP, NSP, nst and nso. The results suggested that the fusion gene was detectable in numerous but not all samples (Table SI-1).

The influence of the presence of the fusion gene on the behavior of the wild-type genes was then examined. The basic expression levels of EWSR1 and further ETS fusion partners and FUS were high throughout and exhibited distinct DE profiles (Fig. S1). The DE profiles support the clinical sample classification and, consequently, the validity of the sample design of the present study. It is worthwhile to observe how the wild-type EWSR1, FUS and ETS genes are influenced by EWSR1-ERG in the CADO cells. Of note, ERG is downregulated in ES2 (EWSR1-FLI1 fusion) but not in ES1 (EWSR1-ERG fusion) samples, but FLI1 expression is not differentially influenced by the fusion gene or the wild-type.

In addition, the dataset was examined for the normalized expression of certain further genes relevant to Ewing sarcoma. Suvà et al (9) mentioned certain markers in the context of Ewing CSC, while Amaral et al (26) analyzed certain markers which are able to segregate Ewing and MSC cells. It appears that in general, these markers also cluster our dataset according to their described properties (Fig. S2A) and reproduce the known origin of our cell types. Another important gene, POLR2A, the RNA polymerase II responsible for transcription and involved in RNA splicing, was reported to interact with EWSR1 (32). The expression profiles of POLR2A together with its interaction partners according to the Search Tool for the Retrieval of Interacting Genes and proteins (STRING) database (https://string-db.org) are presented in Fig. S2B. This gene group has its own coherent profile, which differs from that of the EWSR1-ETS group, but also arranges the samples according to their respective cell types.

Detection of DS by rMATS

To search for alternative splice events between the samples of the experimental design, rMATS was applied (21,33). It was selected as it is still maintained, is able to handle replicated RNA samples (which not all tools support), has a reasonable statistical concept suitable for any technology and preprocessing pipeline, may be limited to known splice variants and is one of the few which is able to detect alternative splicing between two conditions (22). rMATS version 3.2.5 and 4.1.0 with TopHat-aligned BAM files were used. The applied parameter settings were as follows: Paired, read length 35 and anchor length 8. The reference genome for the analysis was Homo sapiens GeneCode version 35 (GRCh38.p13 https://www.gencodegenes.org/human/; gtf format). rMATS was utilized in a fixed mode to detect 5 known and annotated splice events: i) Alternative 3′ splice site (A3SS); ii) alternative 5′ splice site (A5SS); iii) mutually exclusive exons (MXE); iv) retained intron (RI); and v) skipped exon (SE) (Fig. S3).

The five rMATS result tables were provided in two versions. One which evaluates splicing variants with only those reads that span splice junctions (JC) and the other one includes reads that span splice junctions and additionally those reads placed fully on the adjacent, alternatively spliced exon region (reads on target) (JCEC) http://rnaseq-mats.sourceforge.net/). All result sets were generated by the JC approach as well as with the JCEC approach to study the splice variant frequencies across these two approaches.

The raw rMATS results were filtered using a false discovery rate (FDR) of 0.05 as the cutoff. Table I presents the total number of significant genes rMATS generated for the JC results (for JCEC see Table SI-2). All the rMATS result sets represent a DS analysis between two selected cell populations and are therefore inherently differential sets.

Table I DS numbers per splice event from all rMATS analyses.a

Table I

DS numbers per splice event from all rMATS analyses.a

Label	rMATS sample 1	rMATS sample 2	Splice events					Sum (total n=8,564)
			A3SS	A5SS	MXE	RI	SE
C1	NSP 1-3	nst 1-3	46	28	6	69	201	350
C2	NSP 1-3	nso 1-3	37	27	7	36	152	259
B	SP 1-3	NSP 1-3	33	16	10	33	93	185
D	SP 1-3	MSC 1-4	82	59	61	169	428	799
B2	NSP 1-3	MSC 1-4	97	86	61	238	536	1,018
G1	NSP 1-3	ES 1	76	48	147	121	417	809
G2	NSP 1-3	ES 2	172	94	164	216	731	1,377
H1	SP 1-3	ES 1	93	49	155	118	376	791
H2	SP 1-3	ES 2	128	65	155	189	548	1,085
E1	SP 1-3	hFIB 1-4	58	52	61	128	313	612
E2	NSP 1-3	hFIB 1-4	87	69	64	176	441	837
E3	MSC 1-3	hFIB 1-4	58	34	40	53	257	442
Fraction, %			11	7	11	18	52	100

a Because rMATS is a differential approach two sample columns are given below denoting the considered samples. To the right the resulting splice event gene counts are shown. The junction count numbers of the rMATS analyses are given on the basis of a false discovery rate of smaller than 0.05. Percent numbers indicate the fraction of observed splice events. The label column links these experimental approaches to our main study design. rMATS, program for splice variant analysis; MXE, mutually exclusive exons; SE, skipped exon; RI, retained intron; A3SS or A5SS, alternative 3′ or 5′ splice site; DS, differential splicing; SP, side population; NSP, main Ewing population; nst, no stain control; nso, no sort control; ES, Ewing sarcoma; MSC, mesenchymal stem cell; hFIB, human skin fibroblasts.

The preference for the JC approach was based on the following reasons. JC is more stringent and handles only junctions of well-established splice events, while the 'reads on target' dilute the precision of the 'junction read'-defined splice variants, by presumably interfering with further non- annotated splice events or splice sites. This becomes apparent by analyzing the Shashimi plots in more detail. The Bayesian approach of rMATS (33) has high sensitivity and generates numerous results. Looking at the quality of the results, there appears to be a tendency to create false-positives concerning the categorization but also to accept too many low-read count results. In addition, because bulk sequencing data are used in the present study, the averaging of all cells may result in a loss of gene candidates with low expression.

The biological replicates of the samples exhibit a certain variance. To obtain an estimate for the robustness of the results, a variance analysis between SP and NSP replicate combinations was implemented. As indicated in Table SI-3 (yellow and grey boxes), there was a certain variation between the replicate combinations. However, it was not possible to identify a good rationale for outlier exclusion. Therefore, it was decided not to eliminate replicates and proceeded with all of them.

To improve the impact of the rMATS results, the score 'percent spliced-in' (PSI) was applied to certain intermediate results and the final intersections to filter for relevant splice variants. The standard formula inclusion reads divided by the sum of inclusion and exclusion reads was used. The thresholds for sample group one and sample group two are a PSI in the range of 0.1-0.9 for each sample and a difference between the samples of 0.5-0.9. With these stringent criteria, the JC approach yielded 305 PSI-stable DS events, while with the JCEC approach, 787 PSI-stable DS events were obtained. All of the result tables of the rMATS analyses are available in Information S4 (FDR <0.05).

Downstream processing of rMATS results in R

The study design with the different sample combinations and the used acronyms are presented in Fig. 2. The DS approaches of rMATS denoted by connector lines were further refined by set arithmetic procedures applied to generate more specific subsets of DS genes (e.g., B -> BC -> BC′).

Figure 2 Study design. The base cell entities are represented by circles with the respective acronyms. Thick borders around the circles indicate Ewing cell samples with the fusion gene EWSR1-ERG. Colors: CSC SP (orange); main sarcoma population cells (NSP, light orange); MSC (blue). ES1 and ES2 (red) are further Ewing primary cell cultures and hFIB (grey) are an off-target control. The rMATS analyses are denoted by acronyms, e.g., B or C1, and displayed along the connector lines. (A) Analyses presented define the specific differences (BC′) of the Ewing side population with CSC properties (SP) to the main sarcoma cell population (NSP) and to a putative cell of origin (MSC) applying internal controls (nst, nso). (B) Analyses H1,2 and G1,2 define additional SP/NSP differences to two other primary Ewing cell cultures used to form generalized subsets BC′∩ H1,2 and NSP′∩G1,2. (C) Analyses E1-E3 with normal human fibroblasts evaluate the context of SP/NSP differences established in panel A. (D) The differential splice variant sets specific for the three core entities are created by intersections of BC, D and B2 (full circles, black IDs) leading to the partial sets denoted by SP′, NSP′ and MSC′ (blue parts with white labels, bidirectional). The remaining yellow parts (BC′, D′ and B2′, white labels) are indicating differential splicing variants specific solely for the respective comparison. The central intersection (grey) is not considered at all and sparsely populated (JC:4, JCEC:30). SP, side population; CSC, cancer stem cell; NSP, main Ewing population; nso and nst, controls of the main Ewing population; MSC, bone marrow-derived mesenchymal stem cells; ES1 and ES2, further Ewing cell cultures; hFIB, human skin fibroblasts; ∩, intersection operator; ∪, union operator.

Due to the high number of raw data tables, the mathematical software platform R (version 3.4.4 and 4.0.3) was used to perform the evaluation, certain tests and also the Gene Ontology (GO) analysis. The rMATS tables were imported into R. The 'gene symbol' column in conjunction with the 'exon start' column indicates that there are certain genes which possess more than one splice variant. For these cases, a new column joining gene name and physical position was introduced in an additional column named IDnew to be able to distinguish between distinct splice variants. The drawback is that not all analyses may be performed easily with IDnew, such as the GO analysis. To evaluate the suspected error by using only gene names and not the more precise IDnew identifier, all set sizes were analyzed based on gene symbol vs. those based on IDnew. The difference in numbers is going down for the very distinct and filtered sets such as BC′. The difference between gene symbol and IDnew is reasonably low and varying for JC from 3 to 8% (Table SI-4). The GO annotation was run with the gene symbols.

The CADO control groups nso and nst may contain products of potential splice effects, which may occur due to staining and cell sorting procedures. Therefore, splice effects specific for nst and nso were removed. The B set corrected by C1 and C2 sets was termed BC, created by B-(C1∪C2) [∪, union operator]. Table II presents the count of genes prior to and after the correction.

Table II Adjusted SP-NSP differential set BC.a

Table II

Adjusted SP-NSP differential set BC.a

Label	Procedure	Splice events					Sum
		A3SS	A5SS	MXE	RI	SE
B	None	33	16	10	33	93	185
BC	B-(C1∪C2)	15	7	7	20	48	97

a The basis for the subtraction was to build a union of C1 and C2. The union gets removed from B, resulting in a corrected and smaller set called BC. In this manner, it was possible to remove numerous differential splice variants due to staining and sorting stress from the SP-NSP differential approach B. The junction counts are presented here (false discovery rate <0.05). MXE, mutually exclusive exons; SE, skipped exon; RI, retained intron; A3SS or A5SS, alternative 3′ or 5′ splice site; B, intersection of side population and main Ewing population; BC, intersection of B with the nst and nso controls of the main Ewing population; nst, no stain control; nso, no sort control; ∪, union; '-', difference.

The BC, B2 and D differential sets build an overlapping triple, so that certain genes or their respective splice variants are not solely part of only one rMATS differential set but further rMATS differential sets (Fig. 2D). Table III displays the number of alternatively spliced genes specific for the biological entities, specifically discriminating the biological entities (for JCEC see Table SI-5-1). An overview and work- flow of all set theoretical procedures are provided in Fig. S4.

Table III Cell type and differential-specific differential splicing count numbers.a

Table III

Cell type and differential-specific differential splicing count numbers.a

Label	Procedure	Splice events					Sum
		A3SS	A5SS	MXE	RI	SE
SP′	B2-(BC∩D)	3	1	1	5	21	31
NSP′	D -(BC∩B2)	3	2	2	6	11	24
MSC′	BC-(D∩B2)	24	29	21	77	177	328
BC′[SP-NSP]	BC-(B2∪D)	8	4	4	6	16	38
B2′[NSP-MSC]	B2-(BC∪D)	60	53	33	141	319	606
D′[SP-MSC]	D-(BC∪B2)	49	28	36	76	203	392

a Cell type- (SP′, NSP′, MSC′) and the differential-specific (BC′, B2′, D′) splice event counts. MXE, mutually exclusive exons; SE, skipped exon; RI, retained intron; A3SS or A5SS, alternative 3′ or 5′ splice site; SP, side population; MSC, mesenchymal stem cells; NSP, main Ewing population; BC′, B2′, D′, unique sets according to Fig. 2D; ∪, union; ∩, intersection; '-', difference.

Important for a specific downstream analysis is that all subsets are disjunct as far as possible, so that crosstalk between sets is minimized. There is an exclusive specificity for the genes included in the respective subsets. Conceptually, all the interlinked differential situations in the tripartite, as denoted in Fig. 2A, may be discriminated precisely and therefore, this procedure improves the specificity of the single result sets.

As an additional internal control to obtain an estimate for strength of the C1/C2 correction, the difference of the specific set BC′ from the specific set of B′ (both without the proportions of B2 and D) was determined. The result of the C1/C2 correction is reasonable and yielded a 40% smaller set of BC′ compared with B′, and BC′ is a 100% subset of B′. Thus, the expected trends in this set theoretical approach are consistent (Table SI-5-2). Comparison to the further Ewing cell cultures ES1 and ES2 is adding more relevance to the basic DS sets and enables the discrimination of EWSR1-ERG and -FLI1 dependencies of DS genes (Fig. 2B).

As a final control for specificity, the SP, NSP and MSC populations were compared with an unrelated cell type, hFIBs resulting in sets E1/E2/E3 (Fig. 2C). Sets E1-E3 were further intersected with adjacent differential approaches (E1 with D, BC; E2 with BC, B2; E3 with B2, D). In all situations a certain overlap could be observed, although there were also distinct differences (Table SI-6). Considering the huge number of basic biological functions all human cells share, this overlap is expected. On the other hand, the clear differences exemplify these cells as another entity and cross-validate the trends described above.

These primary DS (intersection) results are mechanistically characterized by analyzing the DS genes at the GO level and membership in protein complexes. The GO analysis of DS genes, the association of DS genes with known protein complexes and the sampling procedures to validate the sparse GO term overlap are comprehensively described in Data S1.

Results

Roughly uniform genomic distribution of the five basic splice events

The DS events between the samples of the experimental design (Tables I-III and SI, Fig. 3) were detected by rMATS. The program was utilized to search only for the following known and well-established splice events: i) A3SS; ii) A5SS; iii) MXE; iv) RI; and v) SE. The overall sum of classified splice variants is 8,564 on an FDR level of 0.05.

Figure 3 Samples and conditions. (A) Goal to elucidate whether the characterized SP with CSC features has alternative splicing. (B) Conditions of the experimental design to reveal the difference in alternative splicing of the SP cells to further cell entities. (C) The performed independent differential rMATS experiments. Further details are provided in the Methods section. SP, side population; CSC, cancer stem cell; NSP, main Ewing population; nso and nst, controls of the main Ewing population; MSC, bone marrow derived mesenchymal stem cell; ES1 and ES2, further Ewing cell cultures; hFIB, human skin fibroblasts.

In several DE studies, an uneven whole-genome distribution of DE genes was observed [e.g., (34)]. This is frequently explained by coherently working pathway genes and regulatory neighborhoods or by genomic alterations. To the best of our knowledge, it is unclear whether DS genes in CADO may also exhibit such an uneven genomic pattern.

The number of primary rMATS results was reduced for the selection of hard candidates by applying the score PSI as described in the Methods to 305 (PSI 0.1-0.9 and PSI difference 0.5-0.9). The distribution of all primary or the remaining 305 splice variants on the genome did not reveal any clear preferences for CADO and ES cell samples. Slight accumulations on one or the other chromosomes did not distract from the overall picture of a more uniform distribution. The whole-genome maps of the PSI stable DS variants of the JC and JCEC approach for B, G1, G2, H1 and H2 are provided in supplementary Information S1.

Pronounced segregation of GO terms by splice events

The primary rMATS result genes are bound to certain splice events. The genes of a distinct splice event are bound to their concrete GO terms. All the unique GO terms of each of the five basic splice events used, A3SS, A5SS, MXE, RI and SE, of one sample only exhibit a limited overlap between each other. This is presented in an example for three pairwise splice event comparisons of the CSC SP′ set (Table IV, Fig. 4). The definition of SP′ is provided in the Methods section. Similar albeit slightly less segregated clusters are observed when two different entities are compared, e.g., an SP′ intersection with an MSC′ intersection (Information S2). In addition, the JCEC approach tended toward a somewhat weaker separation of splice variants in comparison to the JC approach. On closer inspection, there were certain exceptions, but the segregation trend was visible throughout numerous comparisons (Information S2).

Figure 4 Splice events are segregating GO terms. Pairwise comparisons of GO terms for three basic splice events occurring in the intersection of SP′ indicate pronounced segregation between different splice events (JC approach). Note the very small overlap for SE vs. RI, SE vs. A3SS and A3SS vs. RI and the low sampling P-values associated with each approach. Columns refer to the different splice events as indicated by acronyms SE, A3SS and RI. 'n' denotes the number of GO terms found for each splice event. Black areas represent GO terms existing in both types, while green areas are unique for each group. White color (space) indicates the absence of the particular GO term. GO, gene ontology; SP′, a unique DS subset of the side population; JC, splice junction counts; RI, retained intron; SE, skipped exon; A3SS, alternative 3′ splice site.

Table IV DS gene pattern for CSC SP and main population.

Table IV

DS gene pattern for CSC SP and main population.

Label	Unique gene names	Overlap
BC′
∩H1	ATP13A3 (194412198), HLCS (36947332), OGT (71536177), POC5 (75707736), TCAIM (44367455), TRA2A (23521706)	EPB41 (29058588), SNHG17 (38422091)
∩H2	KMO (241590011), SNHG1 (62851987), STRADA (63728333), TRDMT1 (17161482)
NSP′
∩G1	POLDIP3 (42601969), RMDN1 (86472348), SMARCB1 (23791755)	ACADVL (7223812), LMO7 (75835151), NOP58 (202278211), TM7SF3 (27003235),
∩G2	AC022826.2 (73825146), SELENBP1 (151369003), ZFAND1 (81721226), ZNF283 (43837052)	TSPAN3 (77055788)

[i] The genes presented may be identified by their name and splice exon start position. BC′∩H1,2: The stable subsets across all three different Ewing samples are defining the generalized DS difference of the CSC side population to the main Ewing populations. NSP′∩G1,2: The common DS genes between the ERG/ERG-based main Ewing populations and the ERG/FLI1 cross-comparison. Round brackets, the genomic position of the splice variant. DS, differential splicing; NSP′, unique subset of main Ewing population; BC′, unique set according to Fig. 2A and D; CSC, cancer stem cell; ∩, intersection.

The validity of the results was determined by several sampling approaches. The global perspective indicates that the probability to sample the combination of GO terms in the joined five splice events of one sample is extremely low (r1; Data S1, Table SII-1). A permutation approach on the unique set of GO terms of the five splice events of one sample indicated that Fig. 4 is stable on a reasonably high significance level (r2; Data S1, Table SII-2). The generalized sampling procedure at the gene set level of the five splice events (r3; Data S1) as a further analysis on how the gene combinations impact the results indicates that the gene sets themselves define the segregation and therefore form slightly different GO term compositions (Table SII-3).

Taken together, these unexpectedly clear results indicate a possible path towards functional subspecialization for different splice variants, and on the other hand, the clear segregation observed supports the robustness of the established GO term system.

Typical marker genes are not DS

There are numerous conventional DE markers defined for the purpose of diagnostics or research in the field of Ewing sarcoma. It is worth seeing if those marker genes may have an overlap with DS genes. Searching the differential rMATS comparisons for well-defined EWSR1-ETS markers, including PAX7 (35) or RUNX2 (36), or MSC markers such as GATA6, TRPC4, FLG or TGM2 (37), or SCA1, CD29, CD44 and CD90 (38), revealed that none of these markers appear to be DS. Consequently, DS of typical markers appears to be less relevant for the regulation of the SP-NSP transitions or any other analyzed comparison and their regulatory differences.

Creation of DS subsets with higher specificity

The major objective of the present study was to determine if the characterized CSC SP of the Ewing sarcoma-derived CADO cell line exhibits differences in its splice variants compared to the NSP from this cell line, as well as to the bone marrow-derived MSCs, the putative cell of origin for Ewing sarcoma.

As a control for possible stress response in the NSP population due to fluorescence-activated cell sorting, its differences from unsorted (both unstained and stained) populations (nst and nso, respectively) were also determined. To obtain an accurate representation of the physiologically relevant difference between SP and NSP, those differences were therefore subtracted from the SP vs. NSP difference (intersection BC). Furthermore, SP vs. NSP differences that also show up in the comparison to the MSC population (intersections D and B2, respectively) were also eliminated, resulting in intersection BC′.

The tripartite comparison SP, NSP and MSC (Fig. 2D) was designed to identify those DS elements of the SP-NSP analysis that are monodirectional (BC′) and therefore, specific for this difference as compared to the bidirectional comparison, cell-specific DS elements (SP′, NSP′), which are shared with further comparisons (D, B2). The SP′ and NSP′ cells exhibited far less bidirectional splice variants than the MSC′ cells.

DS genes characterizing CSC side- and main Ewing sarcoma populations

The search for the subset of generally relevant DS genes within BC′ and NSP′ by comparing with further primary Ewing cell cultures sharing the same fusion gene (EWSR1-ERG) and additionally with a cell line carrying a different ETS fusion gene (EWSR1-FLI1) resulted in a 21 and 16% gene overlap of BC′ to the H1 and H2 intersection sets, respectively. The NSP′ fraction exhibited a considerable overlap of 33% for G1 and 37% for G2. The NSP′∩G1,2 and BC′∩H1,2 sets [∩ intersection operator] are disjunct, supporting the robustness of the study design. The lack of overlap between the BC′∩H1,2 and NSP′∩G1,2 sets suggested that they belong to different regulatory circuits (Tables II and SIII).

The Pearson correlation analysis of the DS genes in Table II according to their basic expression profiles revealed that the DS genes themselves have an unspectacularly uniform gene expression profile across the cell samples and exhibit a weak pattern (Fig. 5A) and form an Ewing cluster. A similar pattern may be observed for the full set of all rMATS DS genes (Fig. 5B).

Figure 5 The DS genes separate the sample groups according to their gene expression status. (A) The DS genes of Table I are forming an Ewing sarcoma cluster in a Pearson correlation of the expression profiles. The heatmap is based on the TMM-normalized RNA quantification table. The Pearson correlation was also used for the 2D hierarchical clustering of rows and columns. The color scale is denoting the correlation values from weak similarity (>0) to strong similarity (<1). The group color scale is combining the biological replicates. (B) The expression profiles of the DS genes do not expose strong differential effects between the genes. The scale is based on log2 expression values. The set index color codes the comparisons identified by their acronyms. 'H overlap' and 'G overlap' describe the common genes between H1,2 and G1,2, respectively. DS, differentially spliced; SP, side population; ES, Ewing sarcoma; NSP, main Ewing population; nso, no sort control; nst, no stain control; MSC, bone marrow derived mesenchymal stem cell; hFIB, human skin fibroblasts.

The defined gene sets had a certain variability in their absolute and splice event-dependent relative gene expression; thus, it is essential to ensure the robustness of the observed differ- ences. The PSI score was applied to establish the relevance of candidates. The remaining candidate genes have an average expression in the observed region of around 3-5 reads per kilobase of transcript, per million mapped reads (RPKM) and above, i.e., well above the threshold of 1 RPKM reported to be a sufficiently sensitive level to detect splice variants (39-41). Applying the stringent PSI criteria as described in the Methods to the DS genes of Table II, several of the questionable candidates described previously were discarded (Table SIV). The particularly stable DS genes of the comparisons BC′∩H1,2 are ATP13A3, EPB41 and for NSP′∩G1,2, they are ACADVL, NOP58 and TSPAN3 (Fig. 6).

Figure 6 PSI stable DS genes of the comparisons BC′∩H1,2 and NSP′∩G1,2. Sashimi plots are used to illustrate the structure of the splice events. The DS genes of BC′∩H1,2 can be seen in (A-C). (A) Indicates gene ATP13A3, (B) gene EPB41 exhibiting an SE splice event and (C) again gene EPB41 with another SE splice event in exactly the same region but not jointly as an MXE event. (D-F) The NSP′∩G1,2 splice events. (D) illustrates the gene TSPAN3, (E) the gene ACADVL and (F) NOP58. All the presented genes passed the PSI criteria. Red and blue color indicate the two rMATS groups which were compared. The numbers on the exon connector lines indicate the number of considered junction reads. Read coverage is indicated by RPKM units on the Y-axis while the X-axis is displaying the genomic position. The two identified splice events are indicated below each pair of graphs in black. A broad diversity of read distributions could be observed in the presented splice events. However, the junction count numbers are well supporting each of the DS events. DS, differentially spliced; MXE, mutually exclusive exons; SE, skipped exon; PSI, percentage spliced in-differential splicing quality criterion; chr, chromosome; RPKM, reads per kilobase of transcript, per million mapped reads.

Analyzing the splice variants of Table SIV, which fail to pass the PSI criteria, further observations may be made. The NSP′∩G1 intersection contains POLDIP3 showing a SE splice event, which may be a misclassified MXE type, well-populated in CADO and not expressed in the Ewing control cell line ES1. RMDN1 is exposing an uncommon A3SS type, which appears to be differential, but a highly expressed region in the case of ES1 was not identified, presumably due to not being part of the annotation. SMARCB1 had two splice variants of A5SS and an RI type. The RI type may be a miss-classified A5SS type. A5SS is clearly present; however, it is expressed in NSP, but not in ES1. This may point to a sporadic weakness of the algorithm erroneously reporting DE as a DS effect (Information S3).

Wild-type EWSR1 was reported by rMATS in several sets (Table SV), but not DS due to low PSI score differences. DS of wild-type ERG was not observed at all and DS of wild-type FLI1, observed one time, missed the set PSI criteria by far. Thus, none of the three molecules have any role in any of the analyzed DS comparisons.

CORUM-based analysis reveals numerous tumor-relevant processes for DS candidates

Well-defined protein complex databases are available, which include i.a. information on gene or protein members of a complex, GO information and links to known diseases associated with complex members. The added value is that numerous direct interaction partners may be identified. Hence, protein complex databases are a rich resource to annotate genes by their functional context. Several characteristics for the DS genes of SP′, NSP′, MSC′ and the sets BC′∩H1,2 and NSP′∩G1,2 from Table II are described below for the JC samples and summarized in Fig. 7.

Figure 7 DS genes are part of important complexes. Certain core mechanisms affected by the DS genes resulting from the analysis are presented. The connector lines denote the sets of the analysis that lead to the described sub-functionality. The mechanisms differentiate the cellular states. DS, differentially spliced; SP′, a unique subset of the side population; EWS, Ewing sarcoma; MSC′, a unique DS subset of the bone marrow derived mesenchymal stem cell; NSP′, a unique DS subset of the main Ewing population; EWS-ERG or EWS-FLI1, cells with one of these fusion genes.

The SP′ set only contains few and small complexes with significant hits, i.a. TERF1 (42) (comprising CDC27, CDC16, ANAPC1, CDC23, ANAPC7, ANAPC5, ANAPC4, ANAPC2, telomere length regulation), DISC1-RHOT1 complex (43), Tankyrin 1-tankyrin 2-TERF1 complex (44) (organization of chromosome structure) and cell division cycle complex (45) (CDC27, CDC16, ANAPC7). Further complexes are involved in DNA binding and translation. Disease information exists for the KEOPS-complex (46) proteins OSGEP, TP53RK and TPRKB, which are involved in cell proliferation, protein translation and DNA-damage-response signaling/ apoptosis (47). Another molecule, POT1, is involved in telomere metabolism (48).

The NSP′ set has far more results than the SP′ set, albeit with moderate enrichment. The functional annotation points to DNA conformation modification, transcriptional repression and DNA modification by acetylation/deacetylation (49). The list of complexes comprises i.a. SIN3-ING1b complex II (50) (i.a. SAP30, SMARCA4, SMARCB1, HDAC1, SMARCC2, HDAC2, SMARCC1, SMARCD1), ALL-1 supercomplex, GPX1-SBP1 complex (51), ubiquitin E3 ligase (BMI1, SPOP, CUL3) and ETS2-SMARCA4-INI1 complex. Disease information exists for SBP1, which is downregulated in several cancer types, GPX1, which is associated with cancer risk and development, MTA1, which is involved in metastasis and cancer formation, and the TWIST complex, which was reported to be required for cancer metastasis in vivo (52).

The MSC′ set is different to the sets mentioned above. The enrichment is higher and includes certain master pathways, such as ribosomal activities (53) and cellular signaling, but no reported link to oncogenes.

The BC′ set (specific for the difference of SP-NSP) resulted in a short, complex list. The GO annotation is primarily associated with posttranslational modification of amino acids, protein targeting, sorting and translocation, vesicular transport, modification by phosphorylation and transcriptional activation. There is certain disease information pointing to cancer, particularly tumor suppressor genes LKB1 (54) and MLL5 described to be co-expressed with OGT (55) and USP7, and the RBL2 complex (56) (DNMT1, E2F4, E2F5, HDAC1, SUV39H1 and RBL2), which is associated with breast cancer.

The D′ set (SP-MSC) is large and similar to the results of MSC′ or NSP′. Here, certain translation and DNA conformation modification-related activities appear. Important in this context are the hits regarding the RNA polymerase II complex and the BAF complex (57) reported to have central roles in Ewing sarcoma (19).

The B2′ set (NSP-MSC) is large and encompasses well-known regulatory functions such as NOP56p-associated pre-rRNA complex, large Drosha complex, COP9 signalosome complex (58), nucleosomal methylation activator complex and DDB2 complex (59). Disease information exists for DGCR8, which is deleted in DiGeorge syndrome, obviously EWSR1, which is involved in Ewing's sarcoma disease and BARD1, which is associated with colorectal adenocarcinoma and attenuates the affinity of BRCA1 and RAD51.

The subsets described in Table II (NSP′∩G1,2, BC′∩H1,2) are small but are associated with noteworthy complexes, which partly already appear in the basic sets above.

BC′∩H1 displays 10 complexes and a GO annotation comprising histone H3-K4 methylation, protein hydroxylation, protein methylation, DNA topological change, cytokinesis, regulation of transcription and the cell cycle. Thus, to a certain extent, the biological functionality is a subset of the sets mentioned above, but more specialized on distinct biological roles. The OGT-TRAK1-TRAK2 complex (60) is associated with hydroxylation, while MLL1/2 complexes (61) stand for histone methylation. The KMT2E-OGT-USP7 complex (62), containing MLL5, emphasizes cell cycle control and has a role in primary cervical adenocarcinomas.

BC′∩H2 exhibits only one complex, the LKB1- STRAD-MO25 complex (54), which does not appear in any of the three other subsets in this section. The GO information points to protein phosphorylation. One of the name-giving genes, LKB1, is known as a tumor suppressor gene. This complex has a role in various sporadic cancers, particularly non-small cell lung cancer.

NSP′∩G1 has 30 different complexes. Of these, numerous complexes have GO terms which point towards DNA conformation modification, and DNA binding and transcription. One of these is again the BAF complex (alias SWI/SNF) (63). Another one is the RNA polymerase II complex in 6 different variants. A further complex is NOP56p-associated pre-rRNA complex, which is involved in ribosome biogenesis.

NSP′∩G2 has only 4 complexes, namely NOP56p-associated pre-rRNA complex, Emerin complex 1 (64), Emerin regula- tory complex and GPX1-SBP1 complex (51). The first three are also part of NSP′∩G1. The Emerin complexes are associated with DNA topological change and transcription. The last one, GPX1-SBP1, is involved in cell redox homeostasis. Disease information for Emerin points to Emery-Dreifuss muscular dystrophy, while GPX1-SBP1 complex is associated with cancer risk and cancer development.

Summarizing the results for the complex-based GO analysis reveals a mechanistic pattern that is strongly connected with the control of cellular core processes involving DNA- and transcriptional processing. Numerous DS genes are linked to cancer processes and allow, at the functional level, the differentiation of the side- and main Ewing sarcoma populations.

Taking all results together, it may be stated that the extracted DS candidates have a good ability to distinguish between the aforementioned CADO subpopulations and furthermore, they are also valid in a broader context of Ewing sarcoma samples of ERG- and FLI1-based fusion genes. A striking feature is the observation that GO terms belonging to a distinct sample exhibit a pronounced segregation between the basic splice events.

Discussion

The main objective of the present study was to extend the available knowledge on Ewing CSCs-characterized in an SP of CADO cells-by studying the appearance, distribution and composition of splice variants on exactly the same dataset that was used in a previous CSC expression study from our group (12). The DS gene sets resulting from the study design developed in the present study are uncovering central modifiers between the tumor biology of CSCs and the main population of Ewing sarcoma cells.

The overall read coverage in the present experiment is sufficient to reliably detect major effects. However, a non-negligible number of spurious observations may have been reduced possibly by using 70× instead of 35× per sample finally mapped reads. While it is difficult to compare this value to other studies with different design, it should be pointed out that the value of 35× of the present study lies well within the range reported elsewhere, 15-80× (65-67). Another aspect concerns the DS detection tools and the quality of their results. While rMATS obtained a fairly good rating in a study by Ding et al (22), the tool has certain weaknesses concerning the selection criteria. A more sophisticated rule set or better implemented Bayesian procedure would improve the specificity for marginal cases.

The observed numbers of DS events for each differential approach expose the differences between the examined cell populations. Less than 200 DS events segregate the CSC SP from the NSP. Much more DS events separate these two populations from the MSC population. However, the SP-MSC number is ~30% smaller than the NSP-MSC number, which points to a somewhat higher similarity of the SP to the MSC population than the main Ewing population. This is consistent with an MSC origin of the CSC SP (2-4), but it cannot be excluded that the SP with its CSC nature, typified here as a DS signature, may also be part of other stem cell lineages.

The relative frequency of the five splice events considered was distinctly different when summed over all splice variants, with abundance values for SE of 52%, A5SS of 7%, A3SS of 11%, RI of 18% and MXE of 11% (cf. Table III). In 2008, Sammeth et al (68) published a calculation of splice event proportions averaged over all three human references (GeneCode, EnsEMBL and RefSeq) and reported the following splice event frequencies: SE, 52%; A5SS, 8%; A3SS, 14%; RI, 6%; and others, 12%. Their values are corresponding very well with those of the present study. In 2019, Lau et al (69) also published statistics with ENCODE RNA-seq data and deter- mined the following frequencies: SE, 65%; A5SS, 7%; A3SS, 11%; RI, 10%; and MXE, 7%, which is in good agreement with the present data except for the considerably higher abundance of RI in the present dataset. However, these other two statistics are based on all known splice variants collected in databases, whereas the present study only identified splice variants that are differential between two conditions. Technical reasons for the higher abundance of RI in the present data may be excluded, as even the subsets reveal similar frequencies. Thus, the higher abundance of RI in the present dataset may be of specific relevance for the (Ewing) cancer situation. RI may result either in functionally different or in dysfunctional protein products [cf. e.g., Monteuuis et al (70)]. However, in each case, the biological networks in which the differentially spliced proteins participate would be affected.

The pure amount of detected splice variants at the gene expression level does not prove that all these splice variants may be translated into proteins [criticism is provided by Tress et al (71), but we tend to agree with Blencowe (72)]. On average, there is a lot of support that the major splice variants will be translated into proteins. Huang et al (73) sophisticatedly demonstrated that at least with modern approaches, all 5 reported gene splice variants of RARRES2 (Chimerin) may also be detected in body fluids at the real protein level. Therefore, it may be assumed that the results of the present study that are at the gene level are mostly retaining their relevance at the protein level.

Another observation from the current rMATS results is that commonly known markers were rarely detected, as exemplified by certain MSC and EWSR1-ETS markers. The reason may be that numerous established molecular markers are developed by DE studies and not DS studies; hence, a DS screening study would rarely expose conventional (DE) markers. DS studies do focus on relative expression ratios instead of absolute expression strength. Therefore, DS studies may be well-suited for developing additional independent marker panels based on alternative splicing patterns.

The comparison of the CADO-CSC SP to the main Ewing sarcoma population (SP-NSP, intersection BC′) together with two further primary Ewing cell cultures enable us to dissect those candidate genes in BC′ that may be generalized to SP-ES1 (both ERG fusion) or SP-ES2 (ERG vs. FLI1 fusion) comparisons. In an analogous way, the Ewing sarcoma NSP-ES1 (both ERG fusion) or NSP-ES2 (ERG vs. FLI1 fusion) comparisons enable us to dissect those candidate genes in the NSP′ that may be generalized to the further two primary Ewing cell cultures in the design of the present study. In each case, the generalization should eliminate those genes specific for CADO and thus narrow down the candidate gene set to those more likely to be of general relevance for Ewing sarcoma progression. Hence, the results are expected to be useful for future biological studies on the pathophysiology of Ewing sarcoma.

The resulting four DS gene groups are representing the DS fingerprint for the CSC and the main Ewing sarcoma population. The biological role of the four PSI criteria passing alternate spliced genes is multifaceted due to their specific selection history. However, their known roles are cancer-associated ones. ATP13A3 is involved in regulatory processes of cation transport across membranes and polyamine transport (74), EPB41 is part of the anaphase (membrane and spindle) processes and may have further roles (75), while ACADVL catalyzes the first step of the mitochondrial beta-oxidation pathway (76), NOP58 is required for the biogenesis of small nucleolar RNAs (U3, U8 and U14) (77) and TSPAN3 (a transmembrane protein) mediates signal transduction events in the course of cell development, growth and motility (78). Taken together, these specific DS genes are useful markers and link to alternate tumor pathways.

Analyzing the association of GO terms with splice events as shown in this study is not done so far as can be seen in Zhang et al (79) for cancer. The presented observation that different splice events are segregating GO terms in a statistically significant way appears to be a novel aspect. This is pointing to a generally assumed systemic feature of a biological network that a set of distinct splice variants stands for a physiological state of the cell. Switching the splice event to another one will change the role of this gene into another functional context. To our knowledge, such a result has not been reported in the literature in this clear appearance by a splice event-based clustering of GO terms. To sum up, splice event patterns form a distinct GO term composition, separating even different cell samples, and may be visualized by a splice event-based GO term clustering.

To date, the alternative splicing of the constitutional exons of EWSR1-ERG or EWSR1-FLI1 themselves has remained largely elusive. While Patócs et al (17) was mainly interested in (breakpoint) variants of the ETS-fusion genes EWSR1-FLI1 and -ERG itself, they were not able to detect any alternative splicing in the constitutional exons of the fusion partners. What they observed was the appearance of up to three different fusion products in one patient sample and the breakpoint variability of the fusion products itself. The present observations confirm known EWSR1-ERG fusions but it was not possible to detect alternate splicing of fusion genes. The wild-type EWSR1, ERG and FLI1 genes appeared to be splice event-free in the present study. However, due to the spurious events reported, it cannot be completely excluded that alternative splicing in a wild-type gene may occur.

In several studies, EWSR1-FLI1 is classified as a typical network hub. Network hubs are molecules that are interacting with numerous molecules in parallel and thus possessing a certain centrality in the biological network structure. For EWSR1-FLI1, the tight interaction with the spliceosome network was reported by Selvanathan et al (80), which led to alternative splicing for CLK1, CASP3, PPFIBP1 and TERT. The former three are also part of the rMATS sets of the present study. The above-mentioned study reported on two further subsets associated with alternative splicing. An expression set with 10 well-validated genes (two in the present B, G2 and H2 sets: HDAC8 and EZH2) and an expression set with nine genes which were assumed to be affected by Selvanathan et al (80). MBNL1 in the latter one also appeared in our analysis and was able to distinguish SP and NSP cells from MSC. The results of Selvanathan et al (80) were obtained from protein studies, array data and partly from RNAseq data, so this may account for certain differences. However, the EWSR1-FLI1 fusion gene they are working on may account for the major differences. The existing overlap between their core genes and the present results suggests that EWSR1-ERG and EWSR1-FLI1 share a certain amount of their molecular regulation.

A subsequent study published by Selvanathan et al (19) in 2019 was centered around the BAF complex and its supporting oncogenic role. One of the well-characterized gene members in the BAF complex is ARID1A (alias BAF250a, SMARCF1). A splice variant of this molecule was not able to be detected in the present study, but ARID1B was observed in SP-MSC as a further BAF member. However, the opener for their study, BAF47 (official SMARCB1), is part of the rMATS intersection sets of the present study. The minor differences from the genes of Selvanathan et al (19) may be explained again by the fact that they focused their study on the EWSR1-FLI1 scenario and not on the EWSR1-ERG variant. However, it appears that the observations of Selvanathan et al (19) again share certain mechanistic aspects with the present analysis. In this context, the work of Boulay et al (81) has a central role by describing that in Ewing sarcoma, the BAF complex is recruited and retargeted by EWSR1 as part of the EWSR1-FLI1 fusion protein, a process which will also likely be valid for EWSR1-ERG.

Selvanathan et al (19) described, in their interaction experiments between the BAF complex and EWSR1-FLI1, further BAF members, such as SMARCA4, SMARCB1, SMARCC1 and SMARCC2. Of these 4 SMARC-members, SMARCA4, SMARCB1 and SMARCC2 (but not SMARCC1) were part of the present analysis. The BAF members detected in the present study may be separated into two categories. ACTB, ACTG1, ACTL6A, SMARCA4, SMARCE1, so six of 13 members of the BAF complex (46%), were obtained in the SP and NSP to ES1/ES2 comparisons. The other group, consisting of ACTB, ACTG1, ARID1B, SMARCA2, SMARCB1 and SMARCC2, was obtained in the SP and NSP to MSC comparisons (Table SVI). Thus, the composition of alternatively spliced BAF members characterizes the differences between CSC side- and main Ewing sarcoma populations as well as from the MSCs. The coherence of the results of Selvanathan et al (19) with the present observations suggests a strong role of the BAF complex for the ETS fusion genes and implies similar work mechanisms for the EWSR1-ERG and -FLI1 fusion genes.

Selvanathan et al (19), somewhat in contrast to Boulay et al (81), speculates on wild-type EWSR1 interactions with the EWSR1-FLI1 fusion gene to explain how wild-type EWSR1, which appears to not bind to the BAF complex directly, may also come into play for the BAF complex-associated mechanisms. The EWSR1 binding by oligomerization as described by Spahn et al (82) may there- fore be a good working hypothesis and the known variability of the molecular splice mechanisms (83,84) supports such a concept. There are further studies describing the interaction of both EWSR1-FLI1 and wild-type EWSR1 with the early spliceosome factors SF1 and SNRPC (32,85) or CCND1 (86), illustrating that the EWSR1-FLI1 fusion gene and wild-type EWSR1 share a certain mechanistic overlap. However, none of the three wild-type genes are part of the DS lists of the present study. Therefore, it is not possible for the present results to provide any support for the mechanisms discussed above.

The above-mentioned DS events of the BAF complex may also be able to tune activation or repression of transcriptional activity by histone modification and chromatin remodeling, as described by Zhu et al (87) for muscle-specific gene expression and Nguyen et al (88,89) for neural stem cell proliferation and embryonic development. The PRC2 complex, which was demonstrated to behave antagonistic to the BAF complex by Nguyen et al (88), was indicated to have three DS genes in the present study (RBBP7, LCOR and MTF2). Overall, these authors concluded that epigenetic modifications appear to primarily trigger differentiation and cell fate decisions by transcriptional control. However, the mentioned histone modifications, despite being able to survive several mitotic events, may also be overruled by other mechanisms such as those of long non-coding RNAs, as indicated by Grote and Herrmann (90). DS events may therefore introduce a further regulatory layer, which may even be a crucial element to trigger a malformed differentiation path from MSCs to Ewing CSCs by modifying the activity of the BAF- and PRC2 complexes. A recent overview focusing on ncBAF (GBAF) but also on further BAF family members comprehensively presented the broad spectrum of action the BAF family of complexes is able to orchestrate (91).

A review by Sveen et al (92) provided an overview of molecular factors associated with 'alternative splicing and cancer hallmarks'. Of these, certain factors were detected in the broader context of the study (34%), but by far not all (Table SVII). VEGFA is described to have an extended exon 8 but was detected to be present in the current study with an optional exon 6. Overall, this review does not really fit the present Ewing sarcoma situation. In contrast to Sveen et al (92), the studies of Selvanathan et al (19,80) show a much better fit with the present results and particularly expose the cancer biology by numerous overlapping candidate genes.

Comparing the DE results with the DS results will foremost raise the conceptional question of what is comparable in these two research situations. DE studies on a gene level are considering expression strength differences between sample groups, while DS studies aim to identify expression ratio differences of gene fragments between sample groups. This may include certain aspects of a DE analysis, but not all and not vice versa. From this perspective, a certain overlap between DE and DS analyses is expected.

The gene overlap between the results of a previous gene expression study by our group (12) and the present results based on the same aligned BAM files ranges from 9 to 30%. Of note, the reported remarkable set of oncogenes, which is differential between SP and NSP cells, is contributing to the largest overlap in the present study (Table SVIII). This is once more underlining the already reported importance of alternative splicing for cancer maintenance in general [e.g., (93)] and particularly for the oncogenes themselves, not as a primary source of cancer progression instead of a secondary booster stabilizing the fusion gene's initiated changes on the splicing machinery itself (19).

Comparing different published DE studies is hampered by a large number of effects ranging from cell culture to experimental technology and may therefore only be suitable for primary effects. However, in contrast to this conservative estimate, remarkable synchronicity may be observed with studies of Hu-Lieskovan et al (94) of 22% and again Selvanathan et al (19) of 20% (Table SVIII).

Finally, it should be noted that splice events are only one aspect of the whole regulatory scenario and genomic alterations as well as epigenetic modifications, and several further regulatory levels together model the outcome of a cellular state.

Concluding on all the different aspects of this discussion, it may be stated that all of the present results build upon established knowledge and are extending the depth of details to a systems biology level.

In conclusion, the comprehensive and systematically structured screening results provided in the present study constitute a valuable data resource for follow-up analyses focusing on specific research questions on Ewing sarcoma and cancer in general. A striking feature is the observation that GO terms and splice event types are building exclusive groups, which may be distinguished and visualized in a GO term-based cluster graph. This feature enables further studies to decode organizational characteristics of the regulation scheme behind splice events.

The analysis revealed a high similarity between EWSR1-ERG and -FLI1 fusion gene variants concerning their participation in the BAF complex. This, together with the defined core sets for CSC-like SP cells, as well the main Ewing sarcoma cells and their DS characteristics, is a good starting point to integrate the understanding of all ETS fusion gene members over time, and thus clarify the mode of action of stem cells in Ewing sarcoma.

A splice variant-aware expression analysis as a standard procedure is expected to provide more comprehensive and fruitful insight into the systems biology of CSCs and their derivatives. It may be helpful to implement this integrative approach manually performed in the present study in a single software tool for detecting both DS and DE. Such a tool would be required to integrate not only different tool chains but also different philosophies, as the DE and DS analyses are answering related but slightly different basic questions, i.e., outcome of regulation vs. regulation itself.

Supplementary Data

Availability of data and materials

All data are freely available from public or institutional data repositories or on request. The transcriptome data are hosted in EBI ArrayExpress (E-MTAB-6067; https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-6067/), public access. CORUM protein complex database (http://mips.helmholtz-muenchen.de/corum/), public access. DAVID database (https://david.ncifcrf.gov/) with public access. Further information (specifically Information S1-S4 and R workspace) is available through the institutional website (http://complex-systems.uni-muenster.de/sinfo.html).

Authors' contributions

Conceptualization, MH and EK; methodology, JM and EK; validation, JM, MH and EK; formal analysis and investigation, JM, MH and EK; writing - original draft preparation, JM, MH and EK; visualization, EK; supervision, EK; project administration, MH and EK; funding acquisition, MH and EK. All authors have read and agreed to the published version of the manuscript and approved the authenticity of the raw data.

Ethics approval and consent to participate

The primary study was performed according to the tenets of the Declaration of Helsinki and approved by the Institutional Review Board (Ethics Committee) of the University of Muenster [multicenter trials E.U.R.O Ewing 99 (EE99 NCT00020566, 12/02/1999) and Ewing 2008 (NCT00987636, 01/10/2009); both registered at https://clinicaltrials.gov]. Written informed consent was obtained from all subjects involved in the study or in the case of children from their respective parents/guardians.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Acknowledgments

The authors would like to thank Ms. Annegret Rosemann, Ms. Dagmar Clemens and Mr. Martin Kiefer from the Department of Pediatric Hematology and Oncology, University Hospital Münster, Germany for their excellent technical work in all types of wet laboratory procedures.

Funding

This study was supported in part by the Federal Ministry of Education and Research (BMBF; grant no. 01GM0869, TranSaRNet; and in part by grant no. 01KT1310, PROVABES). The authors also acknowledge the support from the Open Access Publication Fund of the University of Münster.

Abbreviations:

References