EWS-FLI1 regulates and cooperates with core regulatory circuitry in Ewing sarcoma

Abstract Core regulatory circuitry (CRC)-dependent transcriptional network is critical for developmental tumors in children and adolescents carrying few gene mutations. However, whether and how CRC contributes to transcription regulation in Ewing sarcoma is unknown. Here, we identify and functionally validate a CRC ‘trio’ constituted by three transcription factors (TFs): KLF15, TCF4 and NKX2-2, in Ewing sarcoma cells. Epigenomic analyses demonstrate that EWS-FLI1, the primary fusion driver for this cancer, directly establishes super-enhancers of each of these three TFs to activate their transcription. In turn, KLF15, TCF4 and NKX2-2 co-bind to their own and each other's super-enhancers and promoters, forming an inter-connected auto-regulatory loop. Functionally, CRC factors contribute significantly to cell proliferation of Ewing sarcoma both in vitro and in vivo. Mechanistically, CRC factors exhibit prominent capacity of co-regulating the epigenome in cooperation with EWS-FLI1, occupying 77.2% of promoters and 55.6% of enhancers genome-wide. Downstream, CRC TFs coordinately regulate gene expression networks in Ewing sarcoma, controlling important signaling pathways for cancer, such as lipid metabolism pathway, PI3K/AKT and MAPK signaling pathways. Together, molecular characterization of the oncogenic CRC model advances our understanding of the biology of Ewing sarcoma. Moreover, CRC-downstream genes and signaling pathways may contain potential therapeutic targets for this malignancy.


INTRODUCTION
As a developmental cancer carrying few genetic alterations, Ewing sarcoma is the second most common malignancy of bone and soft tissue predominantly occurring in adolescents and young adults (1). EWS-FLI1, the primary fusion driver, rewires fundamentally the transcriptome of Ewing sarcoma cells (2)(3)(4)(5). Nevertheless, in many cases, EWS-FLI1 requires the cooperation of other transcriptional cofactors, such as WDR5 and CBP/p300, to regulate chromatin modification and gene expression (6)(7)(8). Moreover, EWS-FLI1-targeting transcription factors (TFs), such as MEIS1, NKX2-2, SOX2, and OTX2 (7,(9)(10)(11), are required for the fusion driver to fulfil its oncogenic function in Ewing sarcoma cells. Thus, despite having immense capacity in chromatin regulation, EWS-FLI1 still relies on additional cooperators and mediators to orchestrate gene expression programs in Ewing sarcoma. Hence, comprehensive and unbiased characterization of such partners and mediators is critical for further understanding the biology of Ewing sarcoma.
TFs coordinate cell type-specific transcriptional programs typically through regulating distal cis-regulatory elements, including enhancers and super-enhancers. Although hundreds of TFs are expressed to some extend in any given cell type (12), only a handful appear to be critical for establishing and maintaining cell type-specific gene expression network (13)(14)(15). Notably, this small set of TFs (some-times called master TFs) often form a 'Core Regulatory Circuitry (CRC)' (16)(17)(18), wherein each TF not only selfregulates but also co-regulates each other by directly cobinding to their super-enhancers. Although such CRC transcriptional paradigm has been identified in both normal and neoplastic cell types (19,20), its functionality seems to be particularly critical for developmental tumors in children and young adults who carry few somatic genomic alterations. Indeed, some of the best characterized CRC models are established from neuroblastoma, medulloblastoma, as well as PAX3-FOXO1 + rhabdomyosarcoma, all of which have 'quiet' genomic landscapes compared with adulthood tumors (21)(22)(23). For example, in MYCN-amplified neuroblastoma, a set of master TFs (HAND2, ISL1, PHOX2B, GATA3 and TBX2) assemble a functional CRC to orchestrate the unique gene expression program in this cancer type (22). Moreover, in PAX3-FOXO1-driven rhabdomyosarcoma, the fusion protein exploits super-enhancers to set up a CRC machinery in collaboration with the master TFs (MYOG, MYOD and MYCN); this CRC is addicted by rhabdomyosarcoma cells for survival and proliferation (23). Considering that Ewing sarcoma resembles these developmental cancers in terms of having both few genomic alterations and a single epigenomic driver, we postulated that such oncogenic CRC model might exist in Ewing sarcoma, to cooperate with the transcriptional function of EWS-FLI1. The current study was aimed to identify and characterize this CRC apparatus in Ewing sarcoma, and to elucidate its functional significance in this cancer.

Cell culture
Ewing sarcoma cell lines (A673, SKNMC, EW8 and TC71), human embryonic kidney cells 293T (HEK293T) used in this study were described previously (24,25). Human lung cancer cells A549 and breast cancer cells MCF-7 were also used in this study. Cells were grown in Dulbecco's modified Eagle's medium (DMEM) containing 10% fetal bovine serum (FBS), 100 U/ml penicillin and 100 mg/ml streptomycin, and kept in a humidified incubator at 37 • C with 5% CO 2 . Human Bone Marrow Mesenchymal Stem Cells (hBMSC) were purchased from Guangzhou Saliai stem cell biotechnology, and expanded with hBMSC fetal bovine serum-free medium (Saliai, S-10-011). All cell lines were recently authenticated by short tandem repeat analysis.

Antibodies and reagents
The following antibodies were used in the current study  Supplementary Table S1.

Construction of expression and lentiviral vectors
The pCDH-CMV-Flag-EF1-puro-KLF15 expression vector was amplified based on pCDH-CMV-Flag-EF1-puro vector, and a 3xFLAG-tag was added via PCR. pLKO.1puro or pLKO-Tet-On vectors expressing shRNAs targeting FLI1, KLF15, TCF4 and NKX2-2 were constructed and confirmed by DNA sequencing. The single guide RNAs (sgRNAs) were inserted into lentiviral plasmid encoding dCas9 (pLV hU6-sgRNA hUbC-dCas9-KRAB-T2a-Puro, Addgene, #71236). The sequences of sgRNAs were provided in Supplementary Table S2. To produce viral particles, the recombinant viral vectors and packaging vectors were co-transfected into 293T cells. Supernatants containing viral particles were harvested at 48 h and filtered through a 0.45 M filter after transfection. A673 cells were then infected with the virus in the presence of 10 mg/ml polybrene.

siRNA knockdown
Control scramble and target siRNAs were purchased from Shanghai Genechem Co. Cells were transfected with 100 nM siRNA in OPTIMEM-I (Gibco) for 24-48 h using Lipofectamine RNAiMAX transfection reagent (Invitrogen). Two independent siRNA oligonucleotides were used for each gene.

Chromatin immunoprecipitation (ChIP) and data analysis
ChIP was performed using our previously-described methods with slight modifications (11,26). 2 × 10 7 -2 × 10 8 of A673 cells were harvested and fixed with 1% paraformaldehyde for 10 min at room temperature. The fixation process was terminated by adding 250 mM of glycine. Chromatin solutions including lysis/wash buffer, shearing buffer and dilution buffer were prepared following a standard protocol (18). Samples were washed with PBS and lysed twice with lysis/wash buffer. After filtering through a 29 G needle, samples were harvested by spinning at 13 000 rpm for 10 min at 4 • C. Sample pellets were then resuspended in shearing buffer and sonicated to shear genomic DNA to 300-500 bp. Sonicated sample lyses were subsequently spun at 13 000 rpm for 10 min in 4 • C to remove debris; the supernatants were diluted with dilution buffer. For immunoprecipitation, solubilized chromatin was incubated and rotated with 5 g of anti-KLF15, anti-TCF4, anti-NKX2-2 antibody or IgG control antibody overnight at 4 • C. Antibodychromatin complexes were pulled down by Dynabeads Protein G (Life Technologies) for 4 h at 4 • C. The beads were washed with lysis/wash buffer followed by cold TE buffer. Finally, bound DNAs were eluted by elution buffer (1% SDS, 100 mM NaHCO 3 ) and reverse-crosslinked overnight at 65 • C. DNA molecules were next treated with RNase A and proteinase K. Immunoprecipitated DNAs were extracted with the Min-Elute PCR purification kit (Qiagen), followed by either qPCR analysis, or DNA library preparation and sequencing on the HiSeq 4000 platform (Illumina).
ChIP-Seq was analyzed using established pipelines (11,18). Raw reads were aligned to hg19 reference genome using bowtie 2 aligner (version 2.3.4.3) (27) followed by removal of PCR duplicates with Picard Mark Duplicates (http://broadinstitute.github.io/picard/). ChIP-Seq peaks were called using MACS (Model-Based Analysis of ChIP-Seq, version 2.1.2) (28) with default parameters. Reads were extended in the 5 to 3 direction by the estimated fragment length and normalized at -log 10 of the Poisson P-value of IP file compared to expected background counts by MACS2 bdgcmp command. BigWig files were generated by bedGraphToBigWig tool (https://genome.ucsc.edu/) and visualized in Integrative Genomics Viewer (http://www.broadinstitute.org/igv/home). Motif identification and comparison were performed with HOMER using findMotifsGenome program. H3K27ac ChIP-Seq data generated in four Ewing sarcoma cell lines, three primary tumors and primary mesenchymal stem cells (MSCs) were retrieved from NCBI Gene Expression Omnibus (GSE61953) and processed uniformly. ChIP-Seq data of KLF15, TCF4 and NKX2-2 have been deposited into the GEO under accession number GSE141493, and the matching input file was obtained from our previous published data (GSM2944109).

Cell cycle analysis
Cells were harvested and washed with PBS, followed by fixation with 70% cold ethanol overnight at 4 • C. Cells were washed twice with PBS and stained with propidium iodide. Cell cycle distribution was detected by SONY SA3800 spectral cell analyzer. Data were analyzed by FlowJo 7.6 software (Tree Star).

Apoptosis assay
Cells were double stained with propidium iodide (PI) and Annexin V by using FITC Annexin V apoptosis detection kit (BD Biosciences) according to the manufacturer's instructions. After staining, cells were analyzed using a BD FACSCanto II flow cytometer. Data were analyzed using FlowJo 7.6 software (Tree Star).

RNA extraction and cDNA expression analysis
Total RNA was extracted using RNeasy mini kit (Qiagen). Purified RNA was reverse-transcribed using Maxima H Mi-nus cDNA Synthesis Master Mix (Thermo Fisher). Quantitative real-time qPCR was performed on the AB7300 Detection System (Applied Biosystems, Foster City, CA, USA) using gene-specific primers and Power SYBR Green PCR Master Mix (Applied Biosystems). Expression of each gene was normalized to GAPDH, and quantified using 2-delta (ct) method.

Luciferase reporter assay
Candidate DNA regions (∼500 bp) were PCR amplified and cloned into either pGL3-Promoter firefly luciferase reporter vector or pGL3-Basic luciferase reporter vector (Promega). Constructs were verified by Sanger sequencing. A673 cells were transfected using BioT transfection reagent. A Renilla luciferase control vector was co-transfected as a control for normalization. After 48 h of transfection, the luciferase activity was measured by the Dual-Luciferase Reporter Assay System (Promega).

Xenograft assays in nude mice
All animal experiments were conducted following protocols approved by the Institutional Animal Care and Use Committee (IACUC) of Cedars-Sinai Medical Center. For the study of Orlistat treatment in A673 xenografts, twelve 5-6 weeks old BALB/c-nu female mice (Taconic Bioscience) were subcutaneously inoculated in their dorsal flanks with a suspension of A673 cells (2.0 × 10 6 ). When the tumors grew to 50 mm 3 , mice were injected intraperitoneally with either Orlistat (200 mg/kg/day) or equal volume of vehicle (10% DMSO, 20% cremophor and 70% NaCl). For the study of either KLF15 or TCF4 knockdown, A673 cells were engineered to stably express either doxycycline (DOX)-inducible scrambled shRNA or shRNAs against either KLF15 or TCF4, and were cultured in DMEM supplemented with 10% Tetracycline (Tet)-free FBS (Biological Industries) before implantation. After tumor inoculation, mice were randomized into two groups, and were fed with 2.5 mg/mL doxycycline containing water to turn on the expression of shRNAs. Tumor size and body weight were measured every 2 days. Mice were sacrificed by CO 2 inhalation when the largest tumors were ∼1.5 cm in diameter; tumors were dissected, weighted, and analyzed.

Computational construction of CRC
CRC was computationally constructed using our published methodology (18,29,30) with important modifications. Given the central role of EWS-FLI1 in establishing Nucleic Acids Research, 2020, Vol. 48, No. 20 11437 the transcriptome of Ewing Sarcoma cells, we required that all candidate CRC members have both EWS-FLI1 binding motif and binding peaks in their super-enhancer regions. Therefore, we focused on A673 cell line, a well-characterized Ewing Sarcoma line which had matched H3K27ac and EWS-FLI1 ChIP-Seq results (7,31). Briefly, we first identified 77 super-enhancer-assigned TFs in A673 cell line as shown in Supplementary Table S3. Next, super-enhancer regions assigned to these TFs were extended 500 bp both upstream and downstream, followed by motif scanning with FIMO, for the identification of super-enhancer-assigned TFs. Finally, motif scanning was applied to identify further potential binding sites from other TFs of all auto-regulated TFs in their extended super-enhancer regions which, as mentioned above, must also have EWS-FLI1 binding motif and peaks. Circuitries were then constructed based on all possible fully interconnected auto-regulatory loops.

RNA-Seq and data analysis
RNAs were isolated using miRNeasy RNA isolation kit (Qiagen). We aligned 150 bp paired-end reads to hg19 (UCSC) genome using Kallisto psedo aligner (32). Reads were counted with tximport Bioconductor package (33) and normalized to gene levels by Transcript level abundances (TPM). Differentially expressed genes were identified by DESeq2 package (34) with adjusted P value < 0.05 and absolute log 2 (fold change) > 0.5. RNA-Seq data of siKLF15, siTCF4 and siNKX2-2 are provided in Supplementary Table S4. Gene Ontology (GO) enrichment analysis was performed using ConsensusPathDB (http://cpdb.molgen.mpg. de/) by using all significant downregulated genes as the foreground genes. For those significant downregulated genes related to lipid metabolism or signal transduction, we also found out the associated pathways by performing KEGG pathway enrichment analysis. For GSEA analysis, we used significantly downregulated genes (adjusted P value < 0.05 and log 2(fold change) < -1) in each knockdown of TF as the annotation base and performed enrichment analyses in all expressed genes with mean TPM values > 0.5 in other two TFs. The RNA-Seq data of A673 upon knockdown of each TF have been deposited into the GEO under accession number GSE141493.

Hi-C interactions
The Hi-C data of SKNMC cell line were downloaded from ENCODE database and processed based on the pipeline published by Dixon et al. (31). We extracted the interactions of the loci of KLF15, TCF4 and NKX2-2 at 40 kb revolution and visualized in UCSC genome Browser (https: //genome.ucsc.edu/index.html).

Statistical analysis
Two-tailed Student's t-test was used to evaluate the statistical difference between two groups, while two-way analysis of variance (ANOVA) was applied for multi-group comparisons. All statistical analyses were performed with SPSS 19.0. Log-rank test was used for survival analysis. Differences were considered statistically significant at P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***); n.s., not significant. Diagrams were created by GraphPad Prism 6.

Identification of CRC under the regulation of EWS-FLI1 in Ewing sarcoma
We recently characterized the super-enhancer landscape in Ewing sarcoma and confirmed the essential role of EWS-FLI1 in regulating the epigenome of this cancer (11). As introduced earlier, considering that CRC is particularly important for childhood developmental cancers having few genomic lesions, we postulated that such CRC model might also contribute to regulating Ewing sarcoma transcriptome through either cooperating with or mediating the function of EWS-FLI1. To test this hypothesis, we first sought to identify mathematically master TFs with high interconnectivity through binding to their super-enhancers by motif scanning using our established method (18,29,30). Because of the central role of EWS-FLI1 in establishing the enhancer landscape in Ewing sarcoma, we made significant modifications of the method by requiring that all candidate TFs have both EWS-FLI1 binding motif and binding peaks in their assigned super-enhancers. We initially focused on the A673 cell line, since it is a well-characterized Ewing sarcoma line with available H3K27ac and EWS-FLI1 ChIP-Seq results (7). As a result, a small set (n = 9) of CRC candidates were identified ( Figure 1A), including NKX2-2 and FOS which are known functional cooperators of EWS-FLI1 (7,37). Because CRC factors have high and specific expression in their corresponding cell types, we interrogated the Cancer Cell Line Encyclopedia (CCLE) dataset and noted that, compared with other five CRC candidates (NFATC2, FOS, IRF2, ZBTB7B and MEF2D, Supplementary Figure S1), the expression of 4 candidate TFs (KLF15, NKX2-2, TCF4 and RREB1) showed restricted expression pattern in Ewing sarcoma cell lines (defined as top 5 among all cell types, Figure 1B). However, it should be noted that the specificity of RREB1 expression was overall poor, as most cancer types had comparable mRNA levels.
As a parallel analysis, Pearson correlation coefficient of the mRNA levels of these 9 candidates was determined, considering the co-regulatory relationship between CRC members. Notably, the same set of four factors (KLF15, NKX2-2, TCF4 and RREB1) displayed strong positive correlations with each other ( Figure 1C). In contrast, this correlation pattern was much weaker in other non-Ewing sarcoma bone cancer cell lines ( Figure 1C, lower panel).
Based on these 4 candidates, we next sought to reconstruct functionally the CRC model and to determine their regulatory relationship with EWS-FLI1, by silencing of either EWS-FLI1 or each candidate individually. Knockdown of EWS-FLI1 strongly decreased mRNA expression of KLF15, NKX2-2 and TCF4, but not RREB1 ( Figure  1D). In contrast, the opposite effect was not observed, i.e. knockdown of candidate TFs had no impact on the level of EWS-FLI1, suggesting that the CRC is under the control of EWS-FLI1. Within the four CRC candidates, knockdown of either KLF15, NKX2-2 or TCF4 decreased the expression of the other two ( Figure 1D), but not RREB1. On the other hand, RREB1-silencing did not produce these co-regulatory effects, demonstrating that in Ewing sarcoma cells, KLF15, TCF4 and NKX2-2 together constitute an inter-connected circuitry. Consistent with mRNA expression, silencing EWS-FLI1 decreased protein level of KLF15, NKX2-2 and TCF4 in Ewing sarcoma cell lines ( Figure 1E). Importantly, all of the regulatory effects in CRC were validated at the protein levels with additional siRNAs in both A673 and EW8 cell lines ( Figure 1F).

EWS-FLI1 directly activates three CRC factors
To elucidate the co-regulatory mechanisms between CRC TFs and how they are controlled by EWS-FLI1, ChIP-Seq was performed using specific antibodies against either KLF15, TCF4 or NKX2-2 in A673 cells. Importantly, these three TFs trio-occupied both super-enhancers and promoters of each other, as well as themselves (Figure 2A, Supplementary Figures S2 and S3), forming an interconnected circuitry as predicted by our method. Enrichment of H3K27ac signals at these super-enhancers was observed across all Ewing sarcoma primary tumor samples (n = 3) and cell lines (n = 4). All three distal super-enhancers, but not promoters, were occupied by EWS-FLI1 in both A673 and SKNMC cells (Figure 2A, Supplementary Figure S2), again validating our method. Furthermore, by re-analyzing the publicly-available Hi-C data of SKNMC cells, we confirmed interactions between super-enhancers with promoters in every gene locus, suggesting direct transcriptional coregulation between KLF15, TCF4 and NKX2-2 by binding to each other's super-enhancers. These data also highlight that EWS-FLI1 directly controls the expression of all three CRC members by occupying their super-enhancers ( Figure  2A).
Previous de novo motif analyses elucidate that cisregulatory elements activated by EWS-FLI1 are strongly enriched for GGAA repeats (7). Consistently, prominent enrichment of GGAA repeats were observed in EWS-FLI1 binding sites at all three CRC TFs (Figure 2A, B, Supplementary Figure S2). To further examine whether and how EWS-FLI1 establishes the super-enhancers of CRC factors, we first interrogated chromatin accessibility of normal primary MSCs infected with an EWS-FLI1 expression vector. In the absence of ectopic EWS-FLI1 expression, these genomic regions had negligible ATAC-Seq peaks or H3K27ac signals, indicating little transcriptional activity. Importantly, ectopic expression of EWS-FLI1 converted these regions to super-enhancers with much higher H3K27ac intensity and chromatin accessibility (top 4 tracks in Figure 2B, Supplementary Figure S2). On the other hand, in A673 cells, knockdown of EWS-FLI1 drastically decreased H3K27ac signal in EWS-FLI1-occupied superenhancers (bottom 4 tracks in Figure 2B, Supplementary Figure S2). These results strongly suggest that EWS-FLI1 directly initiates and maintains the super-enhancers of CRC members.
Next, we identified three candidate enhancer constituents (CE1, CE2 and CE3) within the super-enhancer of KLF15 based on the occupancy of TFs ( Figure 2C). Specifically, CE1 and CE3 were trio-occupied by EWS-FLI1, TCF4 and NKX2-2, while CE2 was trio-occupied by KLF15, TCF4 and NKX2-2 ( Figure 2A, C). The spacing of TF binding sites can provide insights into the regulatory organization and cooperativity between TFs (38). We thus performed HOMER motif analysis and observed that the spacing between the binding sites of KLF15, TCF4 and NKX2-2 within CE2 was averaged 23 bp (16-31 bp, Supplementary Figure S2C), supporting the functional cooperation between these TFs. These enhancer elements, as well as a control region, were subsequently cloned into the pGL3promoter luciferase reporter vector which were transfected into A673 cells. Robust reporter activities of CE1 and CE2, but not CE3, were observed ( Figure 2D), consistent with the high H3K27ac signal in CE1 and CE2, but not CE3. We thus consider that CE3 harbors no enhancer activity. To determine the regulation of these enhancers by their occupying TFs, we silenced each TFs, and knockdown of each factor ( Figure 2D) markedly reduced the activity of CE1 and CE2. These results support direct regulation of KLF15 super-enhancer by all CRC members, as well as EWS-FLI1. To validate further the transcriptional functionality of these enhancer elements, we performed CRISPRinterference wherein sgRNAs guide dCas9/KRAB com- plex to suppress targeted cis-regulatory elements (18,39). As shown in Figure 2F, both sgRNAs against CE1 and CE2 significantly decreased H3K27ac modification on these elements (40). Importantly, sgRNAs targeting either CE1 or CE2 markedly reduced the expression of KLF15 ( Figure  2G), as well as the other two CRC TFs (TCF4, NKX2-2) ( Figure 2G), again supporting the interconnected coregulation between these TFs. Moreover, CRC factors, but not the fusion oncoprotein, co-bound the promoter region of KLF15 (Figure 2A). We similarly cloned two promoter regions of KLF15 (PR1 and PR2) into the pGL3-basic luciferase reporter vector, and measured their activities by reporter assays in A673 cells ( Figure 2E). Both of the two promoter regions expectedly showed strong reporter activity, and they were significantly inhibited upon silencing each of the three CRC TFs.

CRC factors cooperate with EWS-FLI1 to orchestrate the transcriptional network of Ewing sarcoma cells
To gain insights into the mechanistic basis of CRC TFs in regulation of the transcriptome of Ewing sarcoma, we analyzed the epigenomic characteristics of their occupancy in A673 cells. To this end, we first annotated the putative promoter (H3K4me3 + /H3K27ac + /H3K4me1 − ) and distal enhancer (H3K4me3 − /H3K27ac + /H3K4me1 + ) regions by available histone marks. Genome-wide peaks of each CRC TFs and EWS-FLI1 were then assigned to these regions. Notably, the occupancy of these TFs was pervasive throughout the genome, with 77.2% of all promoters (8230/10 660) and 55.6% of all enhancers (4496/8087) bound by at least one of the factors ( Figure 3A). Although genomic regions quadruple-bound by all four TFs were rare (0.6% promoters  and 0.4% enhancers), cooperative occupancy (i.e. trio-and dual-occupation) were more common than solo-occupancy in both the promoter and enhancer regions ( Figure 3A, B, Supplementary Figure S3), suggesting strong cooperativity between these factors. Specifically, several important cobinding patterns were observed: (i) validating previous studies (6,7), most EWS-FLI1 peaks were restricted within enhancer regions; (ii) EWS-FLI1 almost always trio-occupied with TCF4/NKX2-2; (iii) most KLF15-binding was in the promoter region, and it either dual-occupied with TCF4 or trio-occupied with TCF4/NKX2-2.
Importantly, regardless of promoter or enhancer elements, trio-binding regions always harbored much higher H3K27ac intensity compared to either solo-or dualbinding regions ( Figure 3C), suggesting that regulatory regions trio-occupied by these TFs may have stronger transcriptional activity. Moreover, tri-occupied regions were more likely to overlap super-enhancers than either dual-or solo-occupied regions ( Figure 3D).
To determine further the possible transcriptional impact of these binding event, matched RNA sequencing (RNA-Seq) data of A673 cells were analyzed. Notably, genes associated with trio-bound promoters exhibited the highest expression levels compared with either solo-or dual-bound promoters ( Figure 3E, upper panel). A similar trend was also observed in enhancer elements, albeit without statistical significance ( Figure 3E, lower panel). Together, these data highlight the cooperative binding as a key characteristic of CRC TFs. Moreover, CRC factors cooperate not only among themselves, but also with EWS-FLI1 in regulating the epigenome of Ewing sarcoma cells.
To establish further the regulatory effect of CRC factors on the transcriptome, we performed RNA-Seq of A673 cells in either the presence or absence of knockdown of each TF. Importantly, gene set enrichment analyses (GSEA) showed that genes decreased following silencing of KLF15 were strongly and significantly enriched in those also downregulated upon depletion of either TCF4 or NKX2-2 (Figure 3F). The same pattern was also observed in RNA-Seq data of either siTCF4 or siNKX2-2 ( Figure 3G, H). Indeed, the downregulated genes in each dataset substantially and significantly overlapped (P < 10 −6 , empirical distribution test); in fact, the 'shared' changes accounted for almost half of all changes in every case (43.8% for TCF4, 51.3% for KLF15, and 47.0% for NKX2-2) ( Figure 3I). These results suggest that CRC factors may act in concert with each other to co-regulate the transcriptome of Ewing sarcoma cells.

CRC factors have strong tumor-promoting function in Ewing sarcoma cells
Considering the prominent roles of CRC in controlling transcriptional network (41,42), we postulated that KLF15, TCF4 and NKX2-2 are required for the viability and proliferation of Ewing sarcoma cells. Indeed, NKX2-2 has been established as a tumor-promoting factor in Ewing sarcoma (43), but the functions of KLF15 and TCF4 remained hitherto unknown in this cancer. We thus performed lossof-function assays and showed that knockdown of either KLF15, TCF4 or NKX2-2 by individual siRNAs ( Figure  4A) markedly inhibited cell proliferation ( Figure 4B, Supplementary Figure S4A) and colony growth ( Figure 4C). The results were verified by doxycycline-inducible expression of multiple independent shRNAs ( Supplementary Figures S4B-E). Fluorescence-activated cell sorting (FACS) analysis showed that depletion of endogenous expression of CRC TFs caused cell cycle arrest ( Supplementary Figures S4F, G). We also detected increased cell apoptosis after knockdown of TCF4 (Supplementary Figure S4H). A previous study showed that NKX2-2 silencing suppressed xenograft growth of Ewing sarcoma (44,45). Here, we tested the dependency of Ewing sarcoma cells on KLF15 and TCF4 in vivo. To this end, we generated cell line models in which either KLF15 or TCF4 expression could be silenced by doxycycline (DOX)-inducible shRNA (A673-TetshKLF15 or A673-TetshTCF4) (Supplementary Figure  S4B, D). Xenograft tumors were then established in nude mice by subcutaneous implantation. Expression of DOXinducible shRNAs against either KLF15 or TCF4 potently inhibited xenograft growth of Ewing sarcoma ( Figure 4D, E). Finally, immunoblotting of xenograft tumors again confirmed the co-regulation of CRC TFs, as shown by the downregulation of each CRC member upon silencing of either KLF15 or TCF4 ( Figure 4F). Taken together, these results demonstrate that as CRC members, KLF15, TCF4 and NKX2-2 may positively co-regulate each other and promote the growth and survival of Ewing sarcoma cells.

CRC factors co-regulate lipid metabolism pathways in Ewing sarcoma
Having established the functional significance of CRC TFs in the proliferation and viability of Ewing sarcoma cells, we next focused on investigating their downstream pathways. Pathway enrichment analysis was performed using downregulated genes, defined as log 2(fold change) < -0.5 and q value < 0.05, upon knockdown of each TF individually. Multiple top enriched terms were shared among each analysis, including 'Signal transduction', 'Metabolism', 'Metabolism of lipids', 'Gene expression (transcription)' and 'RNA Pol-II transcription' (Figure 4G), strongly supporting the functional cooperation between these CRC factors. Among these shared top-ranking pathways, we were particularly interested in the lipid metabolism signaling, since: (i) the biological significance of lipid metabolism in the biology of Ewing sarcoma is unknown; (ii) the regulatory basis of CRC TFs on lipid metabolism is unclear.
To elucidate the mechanisms underlying the regulation of lipid metabolism by CRC TFs, we analyzed in-depth all the enriched genes (n = 242) in the term 'Metabolism of lipids' ( Figure 5A). Again, consistent with the close cooperation between CRC TFs, these genes were strongly shared between individual RNA-Seq data upon knockdown of each CRC factors. These downregulated genes were enriched in key pathways in lipid metabolism, including glycerophospholipid metabolism, sphingolipid metabolism, steroid biosynthesis, biosynthesis of unsaturated fatty-acids, as well as fatty-acid elongation ( Figure 5B). To quantify the functional impact of CRC TFs on these lipid metabolism processes, we performed liquid chromatography-tandem mass spectrometry (LC-MS/MS)-based lipidomics in A673 cells, considering the immense structural complexity of lipid species. Specifically, quantitative lipidomics were performed in the presence and absence of the knockdown of each individual TF. To warrant reproducibility and robustness, 3-4 replicates of each sample were measured, and the data were highly comparably and consistent (Supplementary Figures  S5A-C). In total, we identified an average of 1591 lipid ions in A673 cells, which belonged to 30 lipid classes, suggesting high sensitivity and lipidome coverage of the methodology. The most abundant lipid classes in A673 cells were phosphatidylcholine (PC), phosphatidylethanolamine (PE), diacylglycerol (DG), triacylglycerol (TG), dimethylphosphatidylethanolamine (dMePE) and sphingomyelin (SM). The most abundant fatty acyl chains were oleate (C18:1), palmitate (C16:0), stearate (C18:0), palmitoleic (C16:1), arachidonate (C20:4), docosahexaenoic (C22:6). In agreement with the pathway enrichment analysis, silencing of each TF caused appreciable changes in the lipid landscape of A673 cells ( Figure 5C-E). Specifically, comparing control with knockdown groups, 251 (in siKLF15), 397 (in siTCF4) and 319 (in siNKX2-2) lipid ions were differentially regulated (q value < 0.05, absolute fold change > 2). Although the alterations in lipid species were variable between different knockdown experiments, the majority of which were notably converged to two lipid-associated pathways: glycerophospholipid and sphingolipid pathways (Figure 5F). The numbers on the heatmap show the changed lipids number in each lipid class. The LC-MS/MS-based lipidomics data of siKLF15, siTCF4 and siNKX2-2 in A673 cells are provided in Supplementary Table S5. These changes in lipid classes were highly consistent with and supportive of our pathway analyses of RNA-Seq data, which identified the top two enriched pathways as glycerophospholipid metabolism and sphingolipid metabolism ( Figure  5B). We next examined in detail how lipid metabolism was perturbed by silencing of CRC TFs via integration of RNA-Seq, ChIP-Seq and lipidomic results. This systematic approach identified many important lipid enzymes as direct targets of CRC TFs. For example, rate-limiting enzymes for fatty-acid synthesis and elongation (FASN, ACLY and SCD) were trio-occupied and directly regulated by KLF15, TCF4 and NKX2-2. Moreover, sphingolipid anabolic enzymes (including SPTLC1, DEGS1 and UGCG) and glycerophospholipid anabolic enzyme (LPIN2) were under direct co-regulation by all three CRC TFs ( Figure 5G). Not surprisingly, we also observed solo-or dual-regulation by only 1 or 2 TFs on many other enzymes, such as CERS2 (by KLF15), CERS4 (by TCF4), CERS5 (by NKX2-2) and B4GALT6 (by both KLF15 and NKX2-2). These findings together suggest that CRC factors co-operatively and directly regulate lipid metabolism pathways in Ewing sarcoma cells.

Biological significance of lipid metabolism in Ewing sarcoma cells
Following the characterization of the profound coregulatory effects of CRC factors on lipid metabolic processes, we next sought to explore the biological significance of lipid metabolism in Ewing sarcoma cells. We first tested rate-limiting enzymes for fatty-acid synthesis (FASN, ACLY and SCD, which were trio-regulated by three CRC TFs), since this is the initial process upstream of all lipid metabolism reactions. Verifying our RNA-Seq results, qRT-PCR showed knockdown of each single TF reduced the expression of these central enzymes ( Figure  6A). Immunoblotting confirmed downregulation at the protein levels of FASN, ACLY and SCD ( Figure 6B). In addition, silencing of TCF4 (but not KLF15 or NKX2-2) also reduced the expression of ACC, another key enzyme in fatty-acid synthesis. Some of the trio-binding peaks were shown in Figure 6C using FASN and SCD promotors as examples.
Importantly, silencing of either FASN or SCD reduced cell proliferation and colony growth of both A673 and EW8 cell lines ( Figure 6D, E). Moreover, Orlistat, a specific FASN inhibitor, potently inhibited cell proliferation both in vitro ( Figure 6F) and in vivo ( Figure 6G, Supplementary Figure S5D) where each mouse carried two explants formed by A673 cells, suggesting the biological importance of fattyacid synthesis pathway in Ewing sarcoma cells. Compared with Ewing sarcoma cell lines, human Bone Marrow Mesenchymal Stem Cells (hBMSC) cells were completely resistant to Orlistat ( Figure 6F), suggesting that Orlistat is selective to kill Ewing sarcoma cells. Moreover, although Orlistat inhibited cell viability to varying degrees in other common cancer types, including breast (MCF-7) and lung cancer (A549), they were more resistant to the inhibitor compared with Ewing sarcoma cell lines ( Figure 6F).
In addition to fatty-acid synthesis, we also explored the functional significance of sphingolipid metabolism, another top-enriched lipid pathway in both RNA-Seq and lipidomics ( Figure 5B, G). Among all the sphingolipid anabolic enzymes regulated by CRC, SPTLC1 acts as the ratelimiting enzyme for de novo sphingolipid biosynthesis (46). Indeed, knockdown of each single TF inhibited the expression of SPTLC1 at both transcriptional and protein levels ( Figure 6H, I). Importantly, silencing of SPTLC1 reduced both colony growth and cell proliferation of Ewing sarcoma cells ( Figure 6J-L). Moreover, treatment of Myriocin, a SPTLC1 inhibitor, reduced cell proliferation ( Figure  6M) and colony growth ( Figure 6N) of both A673 and EW8 cells. Similar with the results of Orlistat, hBMSC, MCF-7 and A549 cells were more resistant to Myriocin compared with Ewing Sarcoma cell lines ( Figure 6M).

CRC TFs co-regulate PI3K/AKT and MAPK signaling pathways
In parallel, we also investigated in-depth the term 'Signal transduction' since it was ranked highest in all three RNA-Seq upon knockdown of CRC TFs (1 st in both siKLF15 and siNKX2-2 and 2 nd in siTCF4, Figure 4G). To dissect which specific signaling pathways were coordinately regulated by CRC TFs, KEGG pathway enrichment was next performed using the genes enriched in 'Signal transduction' term (the union set, n = 669) following silencing of each factor ( Figure 7A). Notably, PI3K/AKT and MAPK signaling pathways were identified as the top-ranked pathways ( Figure 7A). Focusing on these two pathways, we found that genes downregulated upon knockdown of each CRC factors were strongly overlapped. Specifically, 32/84 (38%) genes in PI3K/AKT signaling and 23/72 (32%) genes in MAPK signaling were decreased in at least two out of three knockdown groups ( Figure 7B), supportive of the notion that CRC factors cooperatively regulate these pathways.
Further integration with KLF15, TCF4 and NKX2-2 ChIP-Seq data identified that 40/44 (90.9%) commonly downregulated genes (in at least two knockdown groups) of these two cascades were directly occupied by at least one of the CRC TFs ( Figure 7C). Moreover, 23 and 19 of these 40 target genes (57.5% and 47.5%) were dual-and triooccupied by CRC TFs, respectively ( Figure 7C). This result suggests strong cooperation between CRC TFs in the regulation of these two signaling pathways. These direct cotargets included key mediators or effectors of PI3K/AKT and MAPK pathways, such as PDGFRB, PIK3R3, JAK3, VEGFB, RPS6KA5, CCND1 and NFATC1 ( Figure 7C). Indeed, silencing of either KLF15, TCF4 or NKX2-2 inhibited the expression of these key molecules as validated by qRT-PCR ( Figure 7D). Some of the trio-binding peaks were shown in Figure 7E using the RPS6KA5 and PDGFRB promoters as examples. Moreover, immunoblotting assays confirmed reduced phosphorylation levels of several central mediators of MAPK and PI3K/AKT pathways, including P38, ERK, mTOR, P70 and AKT in both A673 and EW8 cells ( Figure 7F). Taken together, these data demonstrate that KLF15, TCF4 and NKX2-2 coordinately promote the transcription of chief components of the PI3K/AKT and MAPK signaling pathways, thereby activating these progrowth and pro-survival signaling cascades in Ewing sarcoma.

DISCUSSION
EWS-FLI1 is a major determinant of genome-wide chromatin states in Ewing sarcoma (2)(3)(4)(5)47,48), by functioning as a pioneer factor (7,9,10,49). Nevertheless, the fusion driver requires other TFs and cofactors to fulfill its oncogenic activities. Considering the essential role of CRC in other developmental cancers driven by single transcriptional regulators (such as NMYC-driven neuroblastoma and PAX3-FOXO1-driven rhabdomyosarcoma) (22,23), we postulated that in Ewing sarcoma, a functional CRC apparatus might cooperate with EWS-FLI1 in the regulation of Ewing sarcoma transcriptome. In order to address this hypothesis, we modified a computational algorithm given the central role of EWS-FLI1, and identified a CRC 'trio' constituted by KLF15, TCF4 and NKX2-2 in Ewing sarcoma cells.
Integrative epigenomics analyses of ChIP-Seq, ATAC-Seq as well as Hi-C together demonstrate that EWS-FLI1, as a pioneer factor, directly establishes the super-enhancers of each of the three CRC TFs to activate their transcription. Subsequently, KLF15, TCF4 and NKX2-2 co-bind to their own and each other's super-enhancers (together with EWS-FLI1 binding) and promoters (without EWS-FLI1 binding), forming an inter-connected auto-regulatory loop ( Figure 7G). However, since knockdown of any individual CRC TFs also reduces other CRC members, it is difficult to fully conclude the functional cooperation between KLF15, NKX2-2 and TCF4. Future investigations including rescue and over-expression assays are required.  In addition to cooperating with EWS-FLI1 to co-amplify their own transcription through a CRC model, KLF15, TCF4 and NKX2-2 also exhibit prominent capability of coregulating the epigenome of Ewing sarcoma cells. In particular, the cooperative occupancy of CRC TFs is pervasive, covering the majority (77.2%) of all putative promoters and over half (55.6%) of all putative enhancers (Figure 3A). Notably, TCF4 and NKX2-2 have higher capacity in occupying the genome than either EWS-FLI1 or KLF15, as their binding regions almost always encompass those of EWS-FLI1 and KLF15. Specifically, TCF4 and NKX2-2 serve as co-binding partners for EWS-FLI1 (in enhancer regions) and for KLF15 (in promoter regions), in addition to their own binding regions. Nevertheless, these three TFs still tend to operate collaboratively, such that their co-bindings events (trio-and dual-) are more common than solo-binding events. Moreover, cooperative binding by more TFs appears to be associated with higher transcriptional activity. Indeed, DNA regulatory elements loaded with more TFs always have higher H3K27ac intensity (Figure 3C), which are also associated with higher expression of downstream genes (only significantly in promoter regions, Figure 3E). These findings are consistent with the proposed function of combinatorial binding of multiple factors in close proximity, which is to overcome with more potency the energetic barrier for nucleosome eviction, thus facilitating activation of cis-regulatory elements.
Not surprisingly, the cooperative binding of CRC TFs results in co-regulation of gene expression program in Ewing sarcoma cells. Indeed, GSEA of RNA-Seq data confirms that downregulated genes upon knockdown of each TF strongly and significantly overlapped with each other. In fact, the shared changes account for almost half of all changes in every RNA-Seq of CRC TFs, strongly suggesting that KLF15, TCF4 and NKX2-2 coordinate to regulate the transcriptome of Ewing sarcoma cells ( Figure 3F-I). Phenotypically, similar to CRC factors in other cancer types, KLF15, TCF4 and NKX2-2 each shows strong progrowth and pro-survival functions in Ewing sarcoma. Indeed, NKX2-2 has been reported to promoting cell proliferation of Ewing sarcoma (43)(44)(45), however, the biological significance of either KLF15 or TCF4 had hitherto been unknown in this cancer.
TCF4 (also known as ITF2) is a basic helix-loop-helix (bHLH) TF that plays a crucial role in the differentiation and specification of central nervous system (CNS) (45,50). Interestingly, depending on different tumor types, TCF4 functions as either an oncogene (diffuse large B-cell lymphoma) (51) or tumor suppressor (medulloblastoma and colon cancer) (52,53). KLF15 is a key regulator of metabolic pathways controlling adipogenesis and gluconeogenesis in both liver and skeletal muscles (54,55). However, its involvement in human malignancies has been poorly understood with mixed findings. For example, KLF15 has shown anti-proliferative activities in gastric and breast cancer cells (56,57); however, KLF15 promotes the proliferation and metastasis of lung adenocarcinoma cells (58). These reports suggest that the cancer-specific functions of KLF15 and TCF4 are highly context-dependent. In our study, consistent with their Ewing sarcoma-promoting functions, all CRC factors are expressed particularly ro-bustly in Ewing sarcoma relative to most of other cancer types ( Figure 1B). In addition to the coordinated function by regulating common genes and pathways, the individual CRC TFs have their own unique biological activities such as cell-cycle progression (Supplementary Figure S4F, G) and apoptosis (Supplementary Figure S4H) by controlling their unique downstream target genes.
To understand the functional contribution of CRC TFs in Ewing sarcoma, we performed pathway enrichment analyses of RNA-Seq upon knockdown of each factors. This approach further substantiates the cooperativity of KLF15, TCF4 and NKX2-2, since multiple top-enriched pathways are shared between three independent loss-of-function experiments. Two of the overlapped high-ranking pathways (lipid metabolism and signal transduction) were prioritized for in-depth investigation given their importance in cancer biology.
Dysregulated lipid metabolism is one of the most important metabolic hallmarks of cancer cells, which is pivotal for synthesis of cellular building blocks and signaling molecules (59)(60)(61)(62). Accumulating evidence suggest that cancer cells depend on altered lipid metabolism for unrestrained growth and survival (62)(63)(64). However, how lipid metabolism is dysregulated and its functional significance in Ewing sarcoma remain poorly understood. In the present study, integrative analyses of LC-MS/MS-based lipidomics, RNA-Seq and ChIP-Seq together highlight that CRC factors converge on regulating key enzymes responsible for the biosynthesis of fatty-acids (FASN, SCD and ACLY), sphingolipids (SPTLC1, DEGS1 and UGCG) as well as glycerophospholipids (LPIN2). These lipid metabolism processes are required for cell survival of Ewing sarcoma, as evidenced by the cytotoxicity of their inhibitors (Orlistat against FASN and Myriocin against SPTLC1). In addition, silencing of either FASN, SCD or SPTLC1 markedly reduced cell proliferation of Ewing sarcoma, further corroborating this conclusion. The biochemical finding that KLF15 regulates lipid synthesis in Ewing sarcoma is not completely surprisingly given the established role of KLF15 in controlling lipid synthesis and fat storage in adipose tissues (65)(66)(67)(68). However, to date, neither TCF4 nor NKX2-2 have been implicated in such metabolic processes in any cell type. These data together highlight CRC-lipid metabolism as a novel protumor cascade in Ewing sarcoma.
As signature oncogenic pathways, PI3K/AKT and MAPK signaling pathways are important for almost every cancer type, including Ewing sarcoma (69)(70)(71). Our investigations ( Figure 7) uncover novel epigenomic mechanisms for the activation of these pathways by CRC factors in Ewing sarcoma. While TCF4 and NKX2-2 have not been involved in these pathways in any cell type, KLF15 is intriguingly reported to inhibit the AKT and MAPK signaling pathways in normal skeletal muscle, cardiomyocytes and kidney (72)(73)(74). It thus appears the regulatory effects of KLF15 on these signaling pathways are cell-type dependent, possibly due to the expression pattern of KLF15 target genes are cell-type specific. This is plausible given that KLF15 almost always cooperates with TCF4 and NKX2-2 in Ewing sarcoma in terms of genomic binding.
In conclusion, we identify and validate an Ewing sarcoma-specific CRC, which is under control of EWS-FLI1. Formed by KLF15, TCF4 and NKX2-2, this CRC apparatus coordinates the gene expression programs, including lipid metabolism, PI3K/AKT and MAPK signaling pathways, in Ewing sarcoma cells. These data advance the understanding of the mechanistic basis of transcriptional dysregulation in Ewing sarcoma, and provide potential novel therapeutic strategies against this malignancy.

DATA AVAILABILITY
ChIP-Seq data of KLF15, TCF4 and NKX2-2, and RNA-Seq data of A673 upon knockdown of each TF have been deposited into the GEO under accession number GSE141493. For ChIP-Seq data, the matching input file was obtained from our previous published data (GSM2944109).

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.