Cooperation of dominant oncogenes with regulatory germline variants shapes clinical outcomes in childhood cancer

Deciphering principles of inter-individual tumor heterogeneity is essential for refinement of personalized anti-cancer therapy. Unlike cancers of adulthood, pediatric malignancies including Ewing sarcoma (EwS) feature a striking paucity of somatic alterations except for pathognomonic driver-mutations that cannot explain overt variations in clinical outcome. Here we demonstrate in the EwS model how cooperation of a dominant oncogene and regulatory variants determine tumor growth, patient survival and drug response. We show that binding of the oncogenic EWSR1-FLI1 fusion transcription factor to a polymorphic enhancer-like DNA element controls expression of the transcription factor MYBL2, whose high expression promotes poor patient outcome via activation of pro-proliferative signatures. Analysis of paired germline and tumor whole-genome sequencing data revealed that regulatory variability at this locus is inherited via the germline. CRISPR-mediated interference with this regulatory element almost abolished MYBL2 transcription, and MYBL2 knockdown decreased cell proliferation, cell survival and tumorigenicity of EwS cells. Combined RNA- and ChIP-seq analyses as well as functional experiments and clinical data identified CCNF, BIRC5 and AURKB as direct MYBL2 targets and critical mediators of its phenotype. In drug-response experiments, high MYBL2 levels sensitized EwS cells for inhibition of its activating cyclin dependent kinase CDK2 in vitro and in vivo, suggesting MYBL2 as a predictive biomarker for targeted anti-CDK2-therapy. Collectively, our findings establish cooperation of somatic mutations and regulatory germline variants as a major determinant of tumor progression and indicate the importance of integrating the regulatory genome in the process of developing new diagnostic and/or therapeutic strategies to fully harness the potential of precision medicine.


INTRODUCTORY PARAGRAPH
Deciphering principles of inter-individual tumor heterogeneity is essential for refinement of personalized anti-cancer therapy. Unlike cancers of adulthood, pediatric malignancies including Ewing sarcoma (EwS) feature a striking paucity of somatic alterations except for pathognomonic driver-mutations that cannot explain overt variations in clinical outcome.
Here we demonstrate in the EwS model how cooperation of a dominant oncogene and regulatory variants determine tumor growth, patient survival and drug response.
We show that binding of the oncogenic EWSR1-FLI1 fusion transcription factor to a polymorphic enhancer-like DNA element controls expression of the transcription factor MYBL2, whose high expression promotes poor patient outcome via activation of proproliferative signatures. Analysis of paired germline and tumor whole-genome sequencing data revealed that regulatory variability at this locus is inherited via the germline. CRISPRmediated interference with this regulatory element almost abolished MYBL2 transcription, and MYBL2 knockdown decreased cell proliferation, cell survival and tumorigenicity of EwS cells. Combined RNA-and ChIP-seq analyses as well as functional experiments and clinical data identified CCNF, BIRC5 and AURKB as direct MYBL2 targets and critical mediators of its phenotype. In drug-response experiments, high MYBL2 levels sensitized EwS cells for inhibition of its activating cyclin dependent kinase CDK2 in vitro and in vivo, suggesting MYBL2 as a predictive biomarker for targeted anti-CDK2-therapy.
Collectively, our findings establish cooperation of somatic mutations and regulatory germline variants as a major determinant of tumor progression and indicate the importance of integrating the regulatory genome in the process of developing new diagnostic and/or therapeutic strategies to fully harness the potential of precision medicine.

MAIN TEXT
The advent of high-throughput 'omics'-technologies in clinical oncology enabled assignment of patients to targeted therapies based on somatic mutations in the protein coding genome 1 .
However, many childhood cancers including Ewing sarcoma (EwS) -a highly aggressive bone-associated cancer -hardly exhibit any recurrent genetic alteration other than pathognomonic and uniformly expressed driver-mutations 2,3 . Yet, these tumors show substantial inter-individual heterogeneity concerning clinical behavior and treatment response, which cannot be solely explained by their few (epi-)genetic alterations [3][4][5] .
However, recent studies in model organisms suggested that the effects of dominant oncogenes may depend on inter-individual variations in the regulatory genome 6 . Thus, we hypothesized that oncogenic cooperation of driver-mutations with specific regulatory germline variants may explain inter-individual diversity in clinical outcomes of oligo-mutated cancers.
We explored this possibility in EwS, which constitutes a genuine model to study such cooperation for several reasons: First, it is characterized by a relatively simple, nearly diploid genome with a single driver-mutation resulting from chromosomal rearrangements fusing the EWSR1 gene to various members of the ETS family of transcription factors (in 85% FLI1) 2,3,7,8 . Second, EWSR1-FLI1 steers ~40% of its target genes by binding DNA at GGAAmicrosatellites, which are thereby converted into potent enhancers [9][10][11] . Third, the enhancer activity of EWSR1-FLI1-bound GGAA-microsatellites strongly depends on the interindividually variable number of consecutive GGAA-repeats 9,10,12 . Together, these characteristics provide an ideal framework to analyze how cooperation of a dominant oncogene (here EWSR1-FLI1) with polymorphic regulatory elements (here GGAA-microsatellites) influences the expression of disease-promoting genes that could explain clinical diversity in oligo-mutated cancers.
To identify such candidate genes with high clinical relevance, we crossed two datasets. The first comprised gene expression microarrays of A673 EwS cells harboring a doxycycline (DOX)-inducible shRNA against EWSR1-FLI1 (A673/TR/shEF1) profiled with/without addition of DOX (Supplementary Table 1 Table 2). We calculated for each gene represented in both datasets the fold change (FC) of its expression after DOX-induced EWSR1-FLI1 knockdown in A673/TR/shEF1 cells and the significance levels for association with overall survival (OS) stratifying patients by expression quintiles of the corresponding gene. This analysis identified MYBL2 (alias B-MYB), encoding a central transcription factor regulating cell proliferation, cell survival and differentiation 13 , as the top EWSR1-FLI1 upregulated gene, whose high expression in primary EwS was significantly associated with poor OS (nominal P=9.6´10 -7 , Bonferroni-adjusted P=0.018) (Fig. 1a,b).
Despite this tight regulation of MYBL2 by EWSR1-FLI1, we noted a marked inter-tumor heterogeneity of MYBL2 mRNA expression in 166 primary EwS (Supplementary Fig. 1e) and in an independent cohort of 208 primary EwS on protein level stained for p-MYBL2 ( Supplementary Fig. 1f). Interestingly, MYBL2 expression did not correlate with minor variations in expression of EWSR1-FLI1 (Supplementary Fig. 1g), suggesting that other factors may cause the inter-individual diversity of MYBL2 transcription in EwS.
In accordance, re-analysis of published chromatin immunoprecipitation followed by nextgeneration sequencing (ChIP-seq) data from A673 and SK-N-MC EwS cells revealed strong signals for EWSR1-FLI1 that mapped to a polymorphic GGAA-microsatellite located ~200 kb telomeric of the transcriptional start site of MYBL2 (Fig. 1d). In both cell lines, this GGAA-microsatellite exhibited EWSR1-FLI1-dependent epigenetic characteristics of an active enhancer indicated by H3K4me1 and H3K27ac marks (Fig. 1d). The EWSR1-FLI1-dependent enhancer activity of this GGAA-microsatellite was confirmed in reporter assays, in which the enhancer activity of cell line-derived haplotypes comprising 6, 10, and 12 consecutive GGAA-repeats was positively correlated with the number of GGAArepeats at this GGAA-microsatellite (Fig. 1e). Applying the haplotype inference and phasing for short tandem repeats (HipSTR) 14 Table 3). Analysis of matched RNA-seq data available for 12 of these 18 EwS primary tumors and one additional tumor sample with WGS and RNAseq data showed that MYBL2 expression was higher in tumors with >15 consecutive GGAArepeats at this microsatellite as compared to tumors with £15 GGAA-repeats (Supplementary Table 3, Supplementary Fig. 1h).
We further validated the EWSR1-FLI1-mediated regulation of MYBL2 in time-course EWSR1-FLI1 ChIP-seq and RNA sequencing (RNA-seq) data generated in A673/TR/shEF1 cells. Removal of DOX after long-term suppression of EWSR1-FLI1 for seven days led to a gradual increase of MYBL2 transcription that correlated with increasing EWSR1-FLI1 recruitment to this GGAA-microsatellite (rPearson=0.816). Strikingly, blockage of protein binding to this MYBL2-associated GGAA-microsatellite by CRISPRi (clustered regularly interspaced short palindromic repeats interference), which promotes an inhibiting chromatin state 15,16 , in highly MYBL2 expressing RDES cells strongly suppressed MYBL2 transcription ( Fig. 1f) and induced a potentially counter-regulatory upregulation of EWSR1-FLI1 ( Supplementary Fig. 1i). Taken together, these findings indicated that MYBL2 is a clinically relevant direct EWSR1-FLI1 target gene, whose expression can be modulated by EWSR1-FLI1 binding to an enhancer-like highly polymorphic GGAA-microsatellite.
To obtain first clues on the functional role of MYBL2 in EwS, we performed gene-set enrichment analysis (GSEA) of MYBL2 co-expressed genes in 166 primary EwS. GSEA revealed that MYBL2 co-expressed genes were significantly enriched in human orthologs of known MYBL2 targets in zebrafish 17 as well as in signatures related to proliferation 18 , cell cycle progression 19 , and sensitization to apoptosis mediated by a CDK-inhibiting protein 20 (min. nominal enrichment score (NES)=2.76, P<0.001) (Fig. 2a, Supplementary Table 4), suggesting that MYBL2 may constitute a key downstream mediator of EWSR1-FLI1-induced, evolutionary conserved proliferation programs.
To test this hypothesis, we performed siRNA-mediated MYBL2 knockdown experiments in three different EwS cell lines (A673, SK-N-MC, and RDES) with moderate to high baseline MYBL2 expression (Supplementary Fig. 2a). Using four different siRNAs ( Supplementary   Fig. 2b,c), we found that knockdown of MYBL2 significantly reduced cell proliferation through blockage of G2/M-progression, which was accompanied by apoptotic cell death ( Fig.   2b-2d). To further explore the function of MYBL2 in EwS growth, we generated two EwS cell lines (A673 and SK-N-MC) with DOX-inducible anti-MYBL2 shRNA expression systems using two different shRNAs (Supplementary Fig. 2d). In both cell lines, DOX-induced MYBL2 silencing significantly reduced clonogenic growth in vitro (Fig. 2e, Supplementary   Fig. 2e) and tumor growth in vivo (Fig. 2f,g, Supplementary Fig. 2f,g) compared to a nontargeting control shRNA. Consistent with our transient knockdown experiments (Fig. 2b-d), we observed an increased number of stalled mitoses, indicating G2/M-blockage, and more apoptotic tumor cells positive for cleaved caspase 3 in xenografts with MYBL2 suppression (Fig. 2h,i). Collectively, these findings indicate that MYBL2 is a critical pro-proliferative downstream effector of EWSR1-FLI1 required for proper G2/M-transition and cell survival.
To identify potential direct MYBL2 targets that could explain its strong pro-proliferative effect, we sequenced RNA of three EwS cell lines with/without siRNA-mediated knockdown of MYBL2 (Fig. 3a). Consistent with our enrichment analyses in primary EwS and functional experiments, GSEA of the identified differentially expressed genes (DEGs) showed that MYBL2 suppression was accompanied by a significant downregulation of the same gene-set comprising human orthologs of zebrafish MYBL2 targets and identical proliferation, cell cycle and sensitization to CDK-inhibitor mediated apoptosis gene signatures (max. NES=-1.89, P<0.001) (Fig. 3b, Supplementary Fig. 3a Fig. 3c). Subsequent ChIP-seq analysis using a specific anti-MYBL2 antibody proved that 66% of these DEGs (50/76) exhibited strong MYBL2 binding at their promoters (Fig. 3c, Supplementary Table 7). Using microarray data of 166 primary EwS, generated with platforms on which 92% (46/50) direct MYBL2 target genes were represented, we correlated the expression levels of these genes with that of Table 8) and these genes' potential association with patients' OS stratifying patients by median expression of the corresponding gene (Supplementary Table 9). Among these genes, CCNF, BIRC5 and AURKB stood out for being highly significantly co-expressed with MYBL2 (Bonferroni-adjusted P<0.05, rPearson≥0.7) (Fig. 3d), and associated with poor OS (P<0.05) (Fig. 3e). To investigate their functional role, we individually knocked down either gene using two specific siRNAs in two different EwS cell lines (Supplementary Fig. 3d) and assessed proliferation and cell viability in vitro.

MYBL2 (Supplementary
Strikingly, knockdown of these genes broadly phenocopied the anti-proliferative and antisurvival effect of MYBL2 silencing (Fig. 3f,g), suggesting that they may constitute important mediators of the pro-proliferative EWSR1-FLI1/MYBL2 transcriptional program.
As there are -to the best of our knowledge -currently no direct MYBL2 inhibitors available, we reasoned that targeting its major upstream cyclin dependent kinase, CDK2, which activates MYBL2 through phosphorylation 13 , may offer a new therapeutic option for EwS patients with high MYBL2 expression. To test this possibility, we treated EwS cells with two different low molecular weight CDK2 inhibitors (CVT-313 and NU6140). While both inhibitors showed strongly reduced growth of A673 EwS cells at the lower micro-molar range ( Fig. 4a), the sensitivity toward both inhibitors was dramatically diminished when MYBL2 was suppressed (Fig. 4a). Such differential effect was not observed in control cells expressing a non-targeting shRNA (Fig. 4a). Treatment of NOD/scid/gamma (NSG) mice with NU6140 significantly reduced growth of EwS xenografts as compared to vehicle (DMSO) (Fig. 4b), and was accompanied by reduced levels of phosphorylated MYBL2 as evidenced by immunohistochemistry (Fig. 4c). However, CDK2-inhibition showed no additional effect on growth of xenografts with silenced MYBL2 expression (Fig. 4b), suggesting that MYBL2 is required for the anti-proliferative effect of CDK2 inhibitors. Consistently, different EwS cell lines with high MYBL2 levels showed higher sensitivity to treatment with NU6140 as compared to a EwS cell line with constitutively low MYBL2 levels ( Supplementary Fig.   4a,b). Since we neither observed significant weight loss (Supplementary Fig. 4c) nor histomorphological changes in inner organs in mice treated for 16 days with up to 40 mg/kg, these results indicated that CDK2-inhibition can safely impair growth of EwS tumors and that MYBL2 may serve as a biomarker to predict their efficacy.
Collectively, our results demonstrate that oncogenic cooperation between a dominant oncogene (here EWSR1-FLI1) and a regulatory germline variant (here a polymorphic enhancer-like GGAA-microsatellite) could be a major source of inter-tumor heterogeneity determining outcome and drug response of an aggressive childhood cancer through modulation of a druggable key downstream player (Fig. 4e). These results also suggest that cooperation between disease-promoting somatic mutations and regulatory germline variants could constitute a general mechanism to explain diversity of disease phenotypes, possibly beyond cancer. In line with this idea, recent reports for neurodegenerative and metabolic diseases showed that the same disease-causing somatic event/mutation can lead to distinct phenotypes depending on (inherited) variations in the regulatory genome 6,21,22 . We anticipate that our findings made in the EwS model are translatable to other malignancies such as ETS-driven leukemia and TMPRSS2-ETS-driven prostate carcinoma 23 , and propose that integration of the regulatory genome in the process of developing new diagnostic and/or therapeutic strategies is necessary to refine and fully exploit 'omics'-based precision medicine.

DNA/RNA extraction, reverse transcription and quantitative real-time PCR (qRT-PCR)
DNA was extracted with the NucleoSpin Tissue kit (Macherey-Nagel), plasmid DNA was extracted from bacteria with the PureYield kit (Promega). RNA was extracted with the NucleoSpin II kit (Macherey-Nagel) and reverse-transcribed using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems). Genomic qRT-PCRs were performed using SYBR green (Applied Biosystems). Oligonucleotides were purchased from MWG Eurofins Genomics. Reactions were run on a Bio-Rad CFX Connect instrument and analyzed using the Bio-Rad CFX Manager 3.1 software. For primer sequences see Supplementary Table 10.

Transient transfection
For siRNA transfection, cells were seeded in a six-well plate at a density of 1.5×10 5 per well in 1.6 ml of growth medium. The cells were transfected with either a negative control non-targeting siRNA (Sigma-Aldrich MISSION siRNA Universal Negative Control #1) or specific siRNAs (25 to 65 nM, depending on the cell line and the siRNA) and HiPerfect (Qiagen). Cells were retransfected 48h after the first transfection and harvested 96h after the first transfection. siRNA sequences are given in Supplementary Table 10.
For plasmid transfection, cells were seeded in a six-well plate at a density of 2×10 5 per well in 1.8 ml of growth medium. Plasmids were transfected with Lipofectamine LTX and Plus Reagent (Invitrogen). The pGL3 vector used for reporter assays has been described before 9 .
Lentivirus production was performed in HEK293T cells. A673 and SK-N-MC EwS cells were infected with the respective lentiviruses and selected with 1.5 µg/ml puromycin (Invivogen).
After single cell cloning, knockdown efficacy of individual clones was assessed by qRT-PCR 48h after addition of DOX (1 µg/ml; Sigma-Aldrich).

DNA constructs and reporter assays
MYBL2-associated microsatellites (with 500 bp 5´ and 3´ flanking regions) from three EwS cell lines were PCR-cloned upstream of the SV40 minimal promoter into the pGL3-luc vector (Promega) 9 . Primer sequences are given in Supplementary Table 10. A673/TR/shEF1 cells (2×10 5 per well) were transfected with the microsatellite-containing pGL3-luc vectors and Renilla pGL3-Rluc vectors (ratio, 100:1) in a six-well plate with 1.8 ml of growth medium.
Transfection media was replaced by media with/without DOX (1 µg/ml) 4h after transfection.
After 72h the cells were lysed and assayed with a dual luciferase assay system (Berthold).
Firefly luciferase activity was normalized to Renilla luciferase activity.

Generation of cell lines with inducible CRISPR interference (CRISPRi)
Due to the lack of functional DNAse, CRISPRi does not cause a knockout of the targeted DNA sequence, but blocks protein binding to it 15,16 . For the reported experiments, a DNAsedead CAS9 (dCAS9) fused to the KRAB effector domain, which promotes an inhibiting chromatin state, is targeted to the genomic region of interest by specific gRNAs to silence the activity of a given enhancer 15,16 . To achieve this, we used a pHAGE TRE dCas9-KRAB vector (Addgene #50917) and a pLKO.1-puro U6 sgRNA BfuAI large stuffer vector (Addgene #52628), the latter one containing either two gRNAs, targeting sequences adjacent to the MYBL2-associated GGAA-microsatellite or a scrambled control (Supplementary Table 10). Lentivirus production was performed in HEK293T cells. RDES EwS cells were infected with the respective lentiviruses and selected with 1 µg/ml puromycin (Invivogen) and 1.5 µg/ml G418 (Invivogen). The cells were induced with DOX for five days, after which MYBL2 and EWSR1-FLI1 levels were measured by qRT-PCR.

Western blot
Protein from A673/TR/shEF1 cells was extracted at d0, d7, d11, d14 and d17 with RIPA and anti-protease cocktail (Roche). Western blots were performed following routine protocols and specific band detection was achieved by the use of rabbit monoclonal anti-FLI1 antibody

Proliferation assays
Cells were seeded in a six-well plate at a density of 1.5×10 5

Analysis of cell cycle and apoptosis
Analysis of cell cycle phases was performed by propidium iodide (PI) (Sigma-Aldrich) staining. Cells were transfected with siRNAs equivalently to the proliferation assays (see above), harvested after 96h (including supernatant), fixed in ethanol (70%) at 4°C, and stained with PI solution (50 µg/ml, with 20 µg/ml RNAse A). Analysis of apoptosis has been performed by combined Annexin V-FITC/PI staining (BD Pharmingen FITC Annexin V Apoptosis Detection Kit II). Cells were transfected with siRNAs equivalently to the proliferation assays (see above) and harvested after 96h (including supernatant). The samples were assayed on an Accuri C6 flow cytometer and analyzed with the Accuri C6 CFlow Plus software. Sample-to-sample normalization and differential expression analyses were performed using the R package DESeq2 (v. 1.18.0) 30 . RNA-seq data were deposited at the Gene Expression Omnibus (GEO; accession code GSE119972).

Chromatin immunoprecipitation and sequencing (ChIP-seq)
DNA-protein cross-linking was done in presence of 1% of paraformaldehyde on 12×10 6 A673 cells for each condition for 10 min. Cell lysis, chromatin shearing, immunoprecipitation and DNA purification was performed with reagents from iDeal ChIP-seq kit for Transcription with option narrow 32 . To normalize, we took the input dataset from the same cell line. PAVIS was used for peak annotation and visualization 33 . ChIP-seq data were deposited at the GEO (accession code GSE119972).

Analysis of published ChIP-seq and DNase-seq data
Publicly available ENCODE SK-N-MC DNase-seq data (GSM736570) 34 and pre-processed A673 and SK-N-MC ChIP-seq data (GSE61944) 35

CDK2 inhibitor assays in vitro
Cells were seeded in a 96-well plate at a density of 5×10 3 per well with/without DOX (1 µg/ml). After 24h of pre-incubation with DOX, CDK2 inhibitors (CVT-313 or NU6140 (Merck)) were added in serially diluted concentrations ranging from 100 nM to 0.001 nM.
Every well contained an equal concentration of 0.5% DMSO. Cells only treated with 0.5% of DMSO served as a control. After 72h of inhibitor treatment, the plates were assayed on a Thermo Fisher Varioskan plate reader after adding Resazurin (20 µg/ml) for 6h.  Survival data were crossed with gene expression microarray data (Affymetrix HG-U133A2.0) generated in A673/TR/shEF1 cells (GSE27524; 53h DOX-treatment) 42 , which were normalized as described above (brainarray CDF, v19).

Construction of tissue microarrays (TMA) and immunohistochemistry (IHC)
The composition and construction of the used TMAs was described previously 44 . All EwS

Evaluation of immunoreactivity
Semi-quantitative evaluation of marker immunostaining was carried out by an independent observer in analogy to scoring of hormone receptor Immune Reactive Score (IRS) ranging from 0-12 as described 45 . The percentage of cells with high nuclear expression was scored in five tiers (score 0=0%, score 1=0-9%, score 2=10-50%, score 3=51-80% and score 4=81-100%). The fraction of cells showing high nuclear expression was scored after examination of 10 high-power fields (40×) of at least one section for each sample. In addition, the intensity of marker immunoreactivity was determined (score 0=none, score 1=low, score 2=intermediate and score 3=strong). The product of both scores defined the final IRS.

ACKNOWLEGEMENTS
We thank Dr. Sarah-Maria Fendt for critical reading of our manuscript, and Anja Heier and Andrea Sendelhofert for excellent technical assistance.

CONFLICT OF INTEREST
We declare no conflicts of interest.