Neuroblastoma, a solid tumor of the peripheral sympathetic nervous system (PSNS) in children, can represent difficult treatment challenges and, as a result, accounts for 15% of all childhood cancer deaths (1, 2). Neuroblastoma is derived from the sympathoadrenal lineage of neural crest cells (NCCs) and typically arises in the adrenal medulla or sympathetic ganglia (3). MYCN gene amplification and overexpression define approximately 25% of neuroblastomas, making it one of the most common high-risk genetic alterations in these tumors (4). The vast majority of neuroblastomas with MYCN gene amplification also harbor large deletions of chromosome band 1p36 (5–8), and a recent study has implicated the gene encoding AT-rich interacting domain–containing protein 1A (ARID1A; also known as BAF250A and SMARCF1) as the gene on chromosome 1p whose loss collaborates with oncogenic MYCN in neuroblastoma pathogenesis (5). This work showed in studies of primary mouse NCCs in vitro that MYCN-driven tumorigenesis selects for NCCs with ARID1A deletions from a pool of NCCs (5).
ARID1A is present exclusively in the canonical BAF complexes within the mammalian SWI/SNF (SWItch/Sucrose Non-Fermentable) (mSWI/SNF) family of chromatin remodeling complexes that modulate chromatin structure and, therefore, gene transcription (9–11). Whole-genome and exome sequencing studies have implicated mSWI/SNF complexes as the most frequently mutated chromatin remodelers in human tumors, with ARID1A being the most highly mutated component among mSWI/SNF subunits (12,13). Inactivating mutations in ARID1A have been identified in a range of tumor types, including neuroblastoma, colon cancer, ovarian clear cell carcinomas, and endometrioid carcinomas (14–19). Furthermore, immunohistochemistry-based studies have demonstrated that most of these mutations are associated with the loss of detectable ARID1A protein expression (20–22).
Sequence alterations of ARID1A have been identified in 6% of neuroblastomas and shown to be associated with early treatment failure and an unfavorable outcome overall (17). ARID1A is deleted on one allele in at least 87% of cases with loss of chromosome 1p, which is almost always deleted in neuroblastomas with MYCN gene amplification and is the most common deletion in high-risk neuroblastomas. The ARID1A gene does not lie within the smallest common region of deletion on 1p, but the vast majority of these abnormalities are very large and include ARID1A within the deleted region (5–8). As a result, at least 30% of all primary high-risk human neuroblastomas have haploinsufficiency for ARID1A (8, 23).
Given the recent studies implicating ARID1A as the critical haploinsufficient tumor suppressor in MYCN-amplified neuroblastoma (5) and the paucity of information regarding the consequences of ARID1A loss in neuroblastoma, we sought to clarify the pathogenic role of this chromatin regulator in our MYCN zebrafish model of high-risk neuroblastoma (24). Here, we show that disruption of either ARID1A homolog, arid1aa or arid1ab, accelerates the onset and increases the penetrance of MYCN-induced neuroblastoma by increasing the proliferation of cells in the sympathoadrenal lineage. ARID1A depletion in human neuroblastoma cells promotes cell invasion and the adrenergic-to-mesenchymal transition by regulating enhancer-mediated gene expression through alteration of the binding sites of both BAF and PBAF (polybromo-containing BAF) complexes (distinct mSWI/SNF assemblies) in a genome-wide manner. These results support the role of ARID1A as a bona fide tumor suppressor in neuroblastoma, whose loss promotes the transition of committed adrenergic neuroblast cells to undifferentiated mesenchymal cells that drive a more aggressive phenotype.
Zebrafish arid1aa or arid1ab deficiency increases the penetrance of MYCN-induced neuroblastoma in vivo
Analysis of gene expression data of human tumors (the R2 database; https://hgserver1.amc.nl/cgi-bin/r2/main.cgi) revealed that low ARID1A expression is strongly associated with lower overall survival probability in neuroblastoma patients and that ARID1A expression levels are inversely correlated with MYCN expression in human primary neuroblastomas (Fig. 1, A and B). To examine the relevance of ARID1A as a tumor suppressor gene in vivo, we used a CRISPR-Cas9–mediated knockout strategy (25) to disrupt its homologs in a zebrafish neuroblastoma model Tg (dβh:EGFP-MYCN) (designated MYCN) (24) and in wild-type (strain AB) fish. Zebrafish have duplicated ARID1A genes, namely, arid1aa and arid1ab. To eliminate potential off-target effects, we designed two unique guide RNAs (gRNAs) targeting the coding sequence of each gene. After injection, randomly selected embryos were evaluated for double-stranded cleavage efficiency of these gRNAs at 1 day postfertilization (dpf). Disruption of the target site for each gRNA was observed in more than 70% of the embryos (table S1). Observation of these mosaic primary injected MYCN fish revealed a substantially accelerated tumor onset and an increased penetrance of neuroblastoma for each gRNA, as compared with uninjected controls (fig. S1). Neuroblastomas arising in these mosaic fish were dissected and subjected to Sanger sequencing to identify gRNA target sites. Mutations in either arid1aa or arid1ab were found in early-onset tumors [5 and 13 weeks postfertilization (wpf)] but not those with late onset (15 and 27 wpf) (table S2), suggesting that the accelerated tumor onset was attributed to CRISPR-Cas–mediated arid1aa or arid1ab gene mutations.
Wild-type fish injected with arid1aa or arid1ab gRNA and Cas9 mRNA were grown to fertility, and stable zebrafish mutant lines arid1aaΔ1 and arid1abΔ19 (arid1aa−/− and arid1ab−/− hereafter) were established by outcrossing (fig. S2, A and B). Each of these mutations included a deletion and/or insertion within a coding region that created a premature stop codon, resulting in a truncation of the Arid1aa or Arid1ab protein before any functional domains, including the DNA binding ARID domain (fig. S2, B to D). Western blotting confirmed the absence of Arid1aa or Arid1ab protein expression in homozygous mutant embryos at 3 dpf (fig. S2E). Zebrafish arid1aa-/-, but not arid1ab-/-, mutants were observed in the adult population. To investigate this result further, we performed quantitative survival studies. While the arid1aa−/− larvae exhibited similar survival rates as wild-type larvae, the arid1ab−/− larvae survival began decreasing at 13 dpf, with no surviving embryos observed beyond 18 dpf (fig. S2F). The arid1ab−/− larvae also displayed a body curvature morphology and a lack of swim bladder indicative of abnormal development (fig. S2G).
To address whether arid1aa or arid1ab deficiency collaborates with MYCN overexpression during neuroblastoma tumorigenesis in the context of stable lines, we incrossed compound heterozygous mutant arid1aa+/−;arid1ab+/− in the background of Tg (dβh:EGFP;dβh:EGFP-MYCN) fish [designated enhanced green fluorescent protein (EGFP); MYCN]. The presence of the dβh:EGFP transgene enhanced fluorescence intensity of the MYCN-fused EGFP, which is more conducive to visualization of tumor development. Zebrafish arid1ab−/− and arid1aa−/−;arid1ab+/− mutants were not detected in the tumor watch population starting at 5 wpf. While fish transgenic for EGFP alone did not develop neuroblastoma regardless of arid1aa and arid1ab genotype (Fig. 1C), both arid1aa or arid1ab deficiency markedly increased the penetrance of MYCN-driven tumors induced in the interrenal gland (IRG; zebrafish counterpart of the human adrenal gland) of EGFP;MYCN fish, with the arid1aa+/−;arid1ab+/−;EGFP;MYCN line showing the highest penetrance: 67.6% in the arid1aa+/−;arid1ab+/−;EGFP;MYCN line versus 28.4% in the control EGFP;MYCN line by 25 wpf (P < 0.0001; Fig. 1C). The mature IRG of EGFP control fish consists mainly of tyrosine hydroxylase (TH)–expressing chromaffin cells (Fig. 1, D to G). By contrast, all MYCN-transformed neuroblastoma cell populations, with or without mutations in the arid1aa or arid1ab gene, included TH-expressing, undifferentiated, small round tumor cells with hyperchromatic nuclei (Fig. 1, H to O), similar to findings in human neuroblastoma (26). Together, these findings demonstrate synergism between arid1aa or arid1ab deficiency and MYCN overexpression in both mosaic and stable zebrafish lines, thus strengthening the case for a tumor suppressor role of arid1aa and arid1ab during neuroblastoma tumorigenesis.
Haploinsufficiency of arid1aa and arid1ab enhances MYCN-induced hyperplasia by increasing cell proliferation in the sympathoadrenal lineage
During normal development, a subset of sympathoadrenal precursor cells migrate ventrally from the neural crest to the anterior-ventral edge of head kidney and differentiate into adrenal chromaffin cells. Sympathoadrenal precursor cells express TH and the panneuronal marker HuC, while chromaffin cells lose HuC expression as they undergo lineage specification, indicating a loss of the neuronal cell identity (24, 27, 28). To clarify the mechanisms responsible for accelerated tumor growth in arid1aa– or arid1ab-deficient MYCN fish, we examined the numbers of TH+HuC+ neuroblasts in sections of the IRG in (i) EGFP, (ii) arid1aa+/−;arid1ab+/−;EGFP, (iii) EGFP;MYCN, and (iv) arid1aa+/−;arid1ab+/−;EGFP;MYCN lines. At 5 wpf, the numbers of TH+HuC+ neuroblasts in the IRG of normal EGFP fish were quite low and the numbers of these PSNS neuroblasts were identical in arid1aa+/−;arid1ab+/−;EGFP fish (Fig. 2A). By contrast, both EGFP;MYCN and arid1aa+/−;arid1ab+/−;EGFP;MYCN fish showed significantly increased numbers of TH+HuC+ neuroblasts in the IRG region, as compared to controls without the MYCN transgene (P = 0.0003; Fig. 2A). The numbers of TH+HuC+ neuroblasts were similar between EGFP;MYCN and arid1aa+/−;arid1ab+/−;EGFP;MYCN fish at 5 wpf (Fig. 2A). At 7 wpf, however, 9 of 10 EGFP;MYCN transgenic fish lacked notable TH+HuC+ neuroblast populations in the IRG region (Fig. 2B) because of a previously reported, developmentally timed apoptotic event (24). By contrast, arid1aa+/−;arid1ab+/−;EGFP;MYCN line continued to display expansion of TH+HuC+ neuroblasts in the IRG at this time point (P = 0.0276; Fig. 2B and fig. S3, G to L), indicating either decreased cell apoptosis or increased cellular proliferation in the sympathoadrenal lineage. arid1aa+/−;arid1ab+/−;EGFP fish showed a slightly increased number of TH+HuC+ neuroblasts in the IRG, as compared to the EGFP control fish at 7 wpf (P = 0.0478; Fig. 2B and fig. S3, A to F), although neuroblastomas were not evident in the lines lacking the MYCN transgene.
To explore this mechanism further, we performed experiments based on activated caspase-3 immunofluorescence and 5-ethynyl-2ʹ-deoxyuridine (EdU) pulse labeling. Apoptotic HuC+ neuroblasts with positive staining for activated caspase-3 were detected in both EGFP;MYCN and arid1aa+/−;arid1ab+/−;EGFP;MYCN fish (P = 1) but not in EGFP and arid1aa+/−;arid1ab+/−;EGFP fish (fig. S3, M to Y). EdU labeling of the proliferating fraction of sympathoadrenal lineage cells in the IRG region of the EGFP;MYCN fish revealed that 4.16% of the TH+HuC+ neuroblasts showed incorporation of EdU after 2 hours of pulse labeling (Fig. 2, C and L to O). By contrast, the percentage of EdU-positive TH+HuC+ neuroblasts was significantly increased to 15.9% in the arid1aa+/−;arid1ab+/−;EGFP;MYCN fish (P = 0.0101; Fig. 2, C and L to S), consistent with the sympathoadrenal hyperplasia observed in this line at 7 wpf (Fig. 2B). No EdU-positive TH+HuC+ neuroblasts were detected in either EGFP or arid1aa+/−;arid1ab+/−;EGFP fish under the same pulse labeling condition (Fig. 2, C to K). Thus, arid1aa+/−;arid1ab+/− haploinsufficiency did not affect the levels of MYCN-induced developmentally timed apoptosis of PSNS neuroblasts but rather increased the numbers of PSNS neuroblasts by increasing the rate of cell proliferation, leading to a higher penetrance of tumors in MYCN fish with arid1aa+/−;arid1ab+/− haploinsufficiency.
ARID1A loss affects enhancer-mediated gene regulation by altering both BAF and PBAF complexes binding in neuroblastoma cells
To assess the tumor-suppressive role of ARID1A in human neuroblastoma cells, we generated two ARID1A homozygous mutant clones of the MYCN-amplified NGP neuroblastoma cell line (designated ARID1A−/− #1 and ARID1A−/− #2) using CRISPR-Cas9 genome editing system with two unique gRNAs (gRNA #1 and #2) targeting human ARID1A (fig. S4). Each of these mutations created a premature stop codon and a truncation of the ARID1A protein product (fig. S4, C and D). An ARID1A heterozygous mutant clone of NGP (designated ARID1A+/−) was also generated using gRNA #2 to assess the consequences of ARID1A haploinsufficiency (fig. S4). A control NGP cell line was generated by transfection with a nontargeting control gRNA. ARID1A protein expression was lower in ARID1A+/− cells and was not detected in either ARID1A−/− #1 or ARID1A−/− #2 cells, when compared to control NGP cells by western blotting (Fig. 3A). Proliferation was unchanged in each of the ARID1A heterozygous and homozygous mutant NGP cells (Fig. 3B). However, transwell invasion and migration assays demonstrated that remarkably more ARID1A−/− #1 and ARID1A−/− #2 NGP cells passed through collagen-coated permeable membranes than did control NGP cells (P < 0.001; Fig. 3, C and D). ARID1A haploinsufficiency also enhanced cell invasion and migration, with an intermediate number of ARID1A+/− NGP cells passing through the membranes compared to control NGP cells (P < 0.01; Fig. 3, C and D). ARID1B is mutually exclusive with ARID1A in BAF complexes and is a specific vulnerability in ARID1A-mutant cancers (29). To test whether these findings extend to neuroblastoma, we first performed an immunoprecipitation experiment with an anti-ARID1B antibody and found that ARID1B-containing BAF complexes were intact in ARID1A knockout NGP cells (Fig. 3E). While knockdown of ARID1B expression with two short hairpin RNAs (shRNAs) did not affect the growth of control NGP cells, it slightly inhibited the growth of ARID1A+/− heterozygous mutant clone and effectively inhibited the growth of both ARID1A−/− #1 and ARID1A−/− #2 homozygous mutant clones (Fig. 3, F and G). Together, these data indicate that ARID1A loss promotes the invasion and migration of human neuroblastoma cells and support the finding of others that ARID1B is synthetic lethal with ARID1A (29), and thus, ARID1B could serve as a therapeutic target in ARID1A-mutant neuroblastomas.
mSWI/SNF family proteins are components of multiple complexes, including the BAF complexes (canonical BAF and recently defined noncanonical BAF complexes) and PBAF complexes (9, 11, 30). Because the genome-wide occupancy of these complexes has not been reported in neuroblastoma, we performed chromatin immunoprecipitation coupled with massively parallel DNA sequencing (ChIP-seq) for BRG1 (shared core subunit in both BAF and PBAF complexes), SS18 (specific in BAF complexes), and ARID2 (BAF200; specific in PBAF complexes) to identify their chromatin localization in control NGP cells. A total of 4346 BAF binding sites defined by overlapping SS18 and BRG1 peaks were identified, while 4017 PBAF binding sites defined by overlapping ARID2 and BRG1 peaks were identified (data file S1 and Materials and Methods). More than half of the PBAF complex binding sites were located at promoter regions, while the BAF complexes preferentially bound to distal intergenic regions (fig. S5A) annotated by ChIPseeker (31). Further analysis showed that about 40% of BAF and PBAF complex binding sites overlapped with each other, and about 60% of the BAF complex target genes were also PBAF complex targets (fig. S5, B and C). Gene ontology (GO) analysis using Enrichr (32) revealed that BAF complex–bound genes in NGP cells were significantly enriched in the categories of nervous system development, cell adhesion, and transcriptional regulation, while PBAF complex–bound genes were mainly enriched for genes involved in transcriptional regulation (fig. S5, D and E).
To investigate the effects of ARID1A loss on the genomic occupancy of BAF and PBAF complexes, we also performed BRG1, SS18, and ARID2 ChIP-seq on ARID1A−/− #1 and ARID1A−/− #2 NGP cells. The BAF or PBAF complex binding sites in control NGP cells were compared to the two homozygous mutant NGP cells, which were computationally merged for downstream analysis (see Materials and Methods for details). A large fraction of BAF binding sites (33.18%) showed weakened SS18 binding in both ARID1A−/− #1 and ARID1A−/− #2 cells compared to control NGP cells, while only a few sites (5.25%) showed strengthened SS18 binding in both ARID1A homozygous mutant cells (Fig. 4A). The genomic occupancy of PBAF complexes was also affected by ARID1A depletion: 11.26% sites gained ARID2-binding capacity, while 4.53% sites lost ARID2 binding in both ARID1A homozygous mutant cells compared to control NGP cells (fig. S6A). Thus, loss of ARID1A not only alters BAF complex targeting of genomic sites but also promotes retargeting of PBAF complexes genome-wide.
To further resolve altered binding activities of the BAF and PBAF complexes and to evaluate the effect of ARID1A deficiency on promoter and enhancer activities, we performed ChIP-seq for H3K4me3 and H3K27ac in all the conditions. Transcription start site (TSS)–proximal/distal BAF and PBAF complex binding sites were classified by distances to the active TSS marked by H3K4me3 (16). We defined TSS-proximal H3K27ac-enriched regions as promoters and TSS-distal H3K27ac-enriched regions as enhancers (see Materials and Methods for details). Thus, we observed widespread changes in H3K27ac levels at enhancers (15.39% up and 27.22% down) but only a few changes at promoters (1.25% up and 2.55% down) upon loss of ARID1A (Fig. 4B). Although a small set of TSS-proximal regions showed altered BAF or PBAF binding, the H3K27ac signal at these sites did not appear to relate to BAF/PBAF binding strength (fig. S7). However, in TSS-distal regions, H3K27ac signal was reduced in both ARID1A−/− #1 and ARID1A−/− #2 homozygous mutant cells at sites with weakened BAF or PBAF binding and was markedly increased in strengthened TSS-distal BAF or PBAF binding sites (Fig. 4, C and D, and fig. S6, B and C). NGP cells with ARID1A haploinsufficiency showed a small effect on BAF but not PBAF binding across TSS-distal regions genome-wide (fig. S8, A and B). Because H3K27ac is indicative of active enhancers and associated with high levels of transcription, we assessed enhancer activities by mRNA expression levels of the nearest genes to BAF/PBAF binding sites. With this approach, changes in the mRNA expression levels of the nearest genes correlated with changes of TSS-distal BAF/PBAF binding upon homozygous loss of ARID1A (Fig. 4E, fig. S6D, and data file S2). GO analysis using Enrichr showed that weakened BAF complex–bound genes upon homozygous loss of ARID1A were selectively enriched in the categories of nervous system development and cell adhesion (fig. S9A), likely indicating these as key pathways involved in ARID1A-mediated tumor suppression. Weakened BAF complex binding sites were selectively enriched for unique transcription factor binding motifs including E-box, aristaless-like homeobox (ALX homeobox), and GATA (fig. S9B). Collectively, these results suggest that ARID1A can both activate and repress enhancer-mediated gene expression in neuroblastoma cells by altering the binding activities of the BAF complexes not only directly but also indirectly because of the retargeting of non-ARID1A containing BAF and PBAF complexes that affect enhancer activity.
ARID1A loss alters the differentiation status of neuroblastoma cells
Neuroblastoma tumors are thought to harbor two unique cell states, adrenergic and mesenchymal, and transcriptional enhancers forming feed-forward autoregulatory or core regulatory circuitry have an important role in programming and reinforcing each cell identity (33–37). We therefore asked whether ARID1A loss would alter the differentiation status of neuroblastoma cells. Our ChIP-seq data showed that paired like homeobox 2B (PHOX2B, an adrenergic-associated transcription factor) gene locus exhibited reduced H3K27ac and BRG1 binding, with fibronectin 1 (FN1, a mesenchymal cell adhesion protein)gene locus showing enrichment with these factors in both ARID1A−/− #1 and ARID1A−/− #2 cells, compared to control NGP cells. Correspondingly, whole-transcriptome mRNA sequencing (RNA-seq) data demonstrated decreased PHOX2B but increased FN1 expression in both ARID1A homozygous mutant cells (Fig. 5A). Western blotting showed that PHOX2B protein expression was completely ablated, while the mesenchymal-associated markers FN1, snail family transcriptional repressor 2 (SNAI2), and vimentin (VIM) were markedly up-regulated in both ARID1A−/− #1 and ARID1A−/− #2 cells compared to control NGP cells (Fig. 5B). RNA-seq analysis of neuroblastoma super-enhancer associated transcription factors (SE-TFs) (37) revealed that 6 of 18 adrenergic SE-TFs had reduced mRNA expression, while 6 of 20 mesenchymal SE-TFs were induced in both ARID1A homozygous mutant cell clones (Fig. 5C). ARID1A+/− heterozygous NGP cells expressed intermediate levels of PHOX2B, FN1, and SNAI2 proteins (fig. S8C). RNA-seq analysis of the same SE-TFs revealed that PHOX2B and GATA2, two major adrenergic SE-TFs, and most of the mesenchymal SE-TFs were also expressed at the intermediate levels in ARID1A+/− heterozygous NGP cells compared to control and ARID1A−/− NGP cells (fig. S8D). Gene set enrichment analysis (GSEA) revealed that a generic gene signature for mesenchymal tumor (38) was significantly enriched in ARID1A homozygous mutant NGP cells (Fig. 5D). Consistent with prior reports of mesenchymal cells being more resistant to conventional chemotherapeutic agents, compared with noradrenergic cells (37), we found that ARID1A homozygous mutant cells were more resistant to cisplatin (Fig. 5E), a standard drug in neuroblastoma therapy. These results provide further evidence that ARID1A loss promotes the transition of NGP cells from an adrenergic to a mesenchymal cell state by regulating enhancer-mediated gene expression.
The ARID1A gene is the most frequently mutated subunit in SWI/SNF complexes in human cancers (12), is deleted on one allele in at least 30% of high-risk neuroblastomas (8, 23), and is mutated in 6% of neuroblastomas overall (17). In this study, we established the first in vivo neuroblastoma model incorporating arid1aa and arid1ab deficiency using zebrafish. Our results demonstrate that loss of just one allele of arid1aa or arid1ab was sufficient to significantly accelerate the onset and increase the penetrance of MYCN-driven neuroblastoma, indicating that loss of even one of the four copies of arid1aa or arid1ab in zebrafish acts to promote the initiation of this tumor. Loss of additional copies of arid1aa or arid1ab further accelerated tumor onset, with the most rapid onset seen in compound heterozygotes of both genes, corresponding to a haploinsufficient tumor suppressor in human cells with only two copies of ARID1A. Our studies emphasize the highly dose-dependent nature of the tumor suppressor role of ARID1A in neuroblastoma tumorigenesis. Although MYCN-induced developmentally timed apoptosis (24, 39) still occurred in the sympathoadrenal lineage cells of arid1aa and arid1ab compound heterozygotes, the increased cell proliferation induced by loss of arid1aa and arid1ab was sufficient to overcome the effects of apoptosis and increase the penetrance of these tumors. Moreover, in studies to further understand the tumor suppressor role of ARID1A in human neuroblastoma, we generated two ARID1A homozygous mutant NGP cell lines, both of which exhibited unchanged cell growth rates (Fig. 3B) but showed increased invasiveness (Fig. 3C) compared to control NGP cells. These results indicate that loss of ARID1A enhances the initiation and migratory properties of neuroblastoma in our model system.
SWI/SNF chromatin remodeling complexes have been shown to be required for the maintenance of lineage-specific enhancers in mouse embryonic fibroblast and rhabdoid tumor cell lines, and SWI/SNF-dependent distal enhancers are essential for controlling expression of genes linked to developmental processes (40). Specifically, SMARCB1 is required for the integrity of SWI/SNF complexes, and its loss alters enhancer activity in rhabdoid tumors (41). Furthermore, ARID1A loss has been reported to impair enhancer-mediated gene regulation in colon cancer (16, 42); similarly, we found that loss of ARID1A reduced genome-wide H3K27ac at enhancer regions in neuroblastoma cells, indicating that they recapitulate this conserved function of SWI/SNF complexes in enhancer regulation. In their study, Mathur and colleagues defined BRG1 and BAF155 cobinding sites as SWI/SNF complexes binding sites, as these factors are shared by both BAF and PBAF complexes (16). By contrast, we considered BAF and PBAF binding sites separately. In wild-type NGP neuroblastoma cells, we found that more than half of all PBAF complex binding sites were located at promoter regions, while BAF complexes preferred to bind at distal intergenic regions (fig. S5A). Our results show that BAF and PBAF complexes colocalize at ~40% of the binding sites and ~60% of the target genes (fig. S5, B and C). Although ARID1A is an exclusive component of the canonical BAF complexes, loss of ARID1A not only impaired BAF complex localization but also affected PBAF complex genomic occupancies of chromatin (Fig. 4 and fig. S6).
SWI/SNF chromatin remodeling complexes are required by leukemia cells for survival and aberrant self-renewal through the maintenance of MYC transcription (43), while the complexes activate p21WAF1/CIP1 through the repression of MYC transcription in differentiation-associated cell cycle arrest (44). These data indicate that the function of SWI/SNF complexes is context dependent. NGP is a MYCN-amplified neuroblastoma cell line, in which the MYCN oncogene is abundantly expressed. Knock out of ARID1A did not significantly alter the expression levels of MYC and MYCN in NGP cells (based on RNA-seq data), indicating that MYC and MYCN are not major downstream effectors of SWI/SNF complexes in the NGP neuroblastoma cell line.
Neuroblastomas include tumor cells exhibiting two different cell states: committed adrenergic cells and undifferentiated mesenchymal or neural crest-like cells, which are more motile and more resistant to chemotherapy. These cells coexist in individual tumor cell populations and exhibit plasticity with the ability to interconvert from one cell state into the other (33, 36, 37). Expression of the mesenchymal transcription factor paired related homeobox 1 (PRRX1) or notch receptor 3 (NOTCH3) intracellular domain can reprogram adrenergic cells toward a mesenchymal state, while PHOX2B, heart and neural crest derivatives expressed 2 (HAND2), and GATA3 work in concert to promote mesenchymal cells to adopt an adrenergic identity (33, 35–37). We found that loss of ARID1A increases the expression of specific mesenchymal markers (FN1, SNAI2, and VIM) and abolishes the expression of an important adrenergic transcription factor dependency gene (PHOX2B), suggesting that ARID1A deficiency promotes the adrenergic-to-mesenchymal transition in neuroblastoma by weakening the influence of BAF complexes that regulate adrenergic enhancer-mediated gene expression (Fig. 5, B and C). Previous studies have shown that increased cell migration and cisplatin resistance in neuroblastoma are accompanied by epithelial-to-mesenchymal transition (37, 45). Thus, the increased invasiveness (Fig. 3C) and cisplatin resistance (Fig. 5E) characteristic of ARID1A mutant NGP cells reflects mesenchymal cell properties acquired upon ARID1A loss. Recent studies to identify and functionally characterize the core regulatory circuitry and signaling feedback loops that reprogram neuroblastoma cells from the adrenergic to the mesenchymal cell state are important (36), as they illuminate pathways that operate in the context of haploinsufficiency for ARID1A to drive the tumor’s metastatic potential and provide a framework for experiments to reveal the most effective therapy.
In conclusion, our work provides compelling in vivo evidence that ARID1A acts as a haploinsufficient tumor suppressor in MYCN-driven neuroblastoma, a prominent and often lethal high-risk tumor in children. Mechanistically, loss of ARID1A has a major impact on the localization of both BAF and PBAF complexes in neuroblastoma cells, with accompanying global effects on transcriptional control. Effects primarily on enhancers of genes are crucial for the determination of cell state and promote adrenergic-to-mesenchymal transition within tumor cell populations of individual tumors. Collectively, these gene expression changes enhance the rate of onset of malignancy and the invasiveness and chemotherapy resistance of the affected tumor cells. Last, our zebrafish model provides an in vivo platform for the identification of small-molecule inhibitors that kill ARID1A-deficient neuroblastoma cells with only minimal toxicity to normal tissues.
MATERIALS AND METHODS
Fish lines and maintenance
Zebrafish were raised and maintained according to the standard procedure. Zebrafish were all AB background strain. All zebrafish studies and maintenance were performed in accord with Dana-Farber Cancer Institute Institutional Animal Care and Use Committee–approved protocol #02-107.
Zebrafish lines Tg(dβh:EGFP) and Tg (dβh:EGFP-MYCN) (24) were described previously and were designated EGFP and MYCN in the text, respectively. The zebrafish arid1aa and arid1ab mutant lines were generated by the CRISPR-Cas genome editing system using the previously described protocol (25). The plasmid constructs pDR274 (#42250) and pMLM3613 (#42251) were purchased from Addgene. Oligonucleotides listed in table S3 were annealed and cloned into the pDR274 vector to generate gRNAs targeting arid1aa and arid1ab. All gRNAs were transcribed in vitro from Dra I linearized pDR274-gRNA plasmid DNA using the MAXIscript T7 Transcription Kit (Ambion, #AM1312). Cas9 mRNA was transcribed in vitro from Pme I linearized pMLM3613 DNA, using the mMESSAGE mMACHINE SP6 Transcription Kit (Ambion, #AM1340). Each embryo was injected with 1 nl of solution containing gene-specific gRNA (25 ng/μl) and Cas9 mRNA (600 ng/μl) at the one-cell stage. Mosaic F0 fish with germline mutation were identified, and stable mutant lines for arid1aa E1-2 gRNA and arid1ab E1-1 gRNA were established by outcrossing. To genotype the arid1aa or arid1ab mutant line, DNA fragment containing the mutated site was amplified with genomic DNA extracted from fin clips serving as template using the primers listed in table S4. The DNA product was either sent for Sanger sequencing with the primers listed in table S4 or digested with the specific restriction enzymes listed in table S4. The DNA product amplified from wild-type allele will be cut into two fragments, while the one amplified from mutant allele will be resistant to digestion.
Mosaic primary injected or stable fish lines were observed biweekly for evidence of EGFP fluorescent tumors starting at 5 wpf. Fish with tumors were separated, genotyped, and euthanized after anesthesia and characterized by histopathological hematoxylin and eosin (H&E) staining and molecular analysis. The cumulative frequency of neuroblastoma development was analyzed by the Kaplan-Meier method with the log-rank test used to determine statistical significance.
Protein extraction, western blot analysis, and immunoprecipitation
Protein was extracted from cell lines or zebrafish samples using radioimmunoprecipitation assay buffer. Western blotting was performed as described previously (46). Protein samples were separated by an 8% polyacrylamide gel and transferred to a polyvinylidene difluoride membrane (Thermo Fisher Scientific, #PI88518). After incubation with antibodies, the results were visualized on autoradiography films with SuperSignal West Pico Chemiluminescent Substrate (Pierce, #PI34080) or SuperSignal West Dura Extended Duration Substrate (Pierce, #34075). Zebrafish Arid1ab was detected with an antibody produced by Josman LLC. that recognizes a region between residue 1770 and 1845 of zebrafish Arid1ab. Anti-ARID1A rabbit polyclonal antibody was purchased from Novus Biologicals (#NB100-55334) and used to detect both human ARID1A and zebrafish Arid1aa. Other antibodies and their sources were anti–α-tubulin mouse monoclonal (B-5-1-2) antibody (Sigma-Aldrich, #T6074), anti-PHOX2B mouse monoclonal (C-3) antibody (Santa Cruz Biotechnology, #SC-376993), anti-SNAI2 rabbit monoclonal (C19G7) antibody and anti-VIM rabbit monoclonal (D21H3) antibody (Cell Signaling Technology, #9585S and #5741S), and anti-FN1 sheep polyclonal antibody (R&D Systems, #AF1918).
Immunoprecipitation was performed as described previously (12). Cells were first lysed in buffer A [10 mM Hepes (pH 7.6), 25 mM KCl, 1 mM EDTA, and 10% glycerol, supplemented with fresh 1 mM dithiothreitol (DTT) and complete protease inhibitors (Roche, #04693132001)] to eliminate cytoplasmic protein. Nuclei were resuspended in buffer C [10 mM Hepes (pH 7.6), 3 mM Mg2Cl, 100 mM KCl, 0.1 mM EDTA, and 10% glycerol, supplemented with fresh 1 mM DTT and complete protease inhibitors] and lysed for 30 min under rotation after the addition of ammonium sulfate to a final concentration of 0.3 M. Soluble nuclear proteins were separated from insoluble chromatin fraction by ultracentrifugation (100,000g) followed by precipitation with ammonium sulfate (0.3 mg/ml). After ultracentrifugation (100,000g), protein precipitate was resuspended in immunoprecipitation buffer [150 mM NaCl, 50 mM tris-HCl (pH 8.0), 1% NP-40, and 0.5% deoxycholate, supplemented with fresh 1 mM DTT and complete protease inhibitors] for immunoprecipitation. The anti-ARID1B mouse monoclonal (KMN1) antibody (Santa Cruz Biotechnology, #sc-32762 X) coupled with Dynabeads M-270 Epoxy (Life Technologies, #14301) was used for immunoprecipitation. Proteins were detected by the anti-ARID1A rabbit polyclonal antibody (Novus Biologicals, #NB100-55334), anti-ARID1B mouse monoclonal antibody (Abcam, #ab57461), anti-BRG1 mouse monoclonal (EPNCIR111A) antibody (Abcam, #ab110641), anti-BAF155 rabbit monoclonal (D7F8S) antibody (Cell Signaling Technology, #11956), anti-SS18 rabbit monoclonal (D6I4Z) antibody (Cell Signaling Technology, #21792), anti-BAF47 mouse monoclonal (F-4) antibody (Santa Cruz Biotechnology, #sc-166164), anti-plant homeodomain finger protein 10 (PHF10) rabbit polyclonal antibody (GeneTex, #GTX116314), and anti–β-actin rabbit polyclonal antibody (Abcam, #ab8227).
EdU pulse labeling, cryosectioning, paraffin sectioning, and immunostaining
For EdU pulse labeling experiment in zebrafish, Edu from the Click-iT Alexa Fluor 647 imaging kit (Life Technologies, #C10340) was diluted to 2.5 mg/ml for retro-orbital injection. Each euthanized fish examined at 5 wpf was injected with 1 μl of EdU solution. Fish were fixed 2 hours after injection for cryosectioning; EdU detection was performed according to the manufacturer’s protocol.
Zebrafish cryosectioning and immunofluorescence staining were performed as previously described (39). Fish at indicated stage were fixed in 4% paraformaldehyde at 4°C overnight, followed by washing three times with phosphate-buffered saline (PBS). The samples were embedded in 1.5% agarose, melted in 30% sucrose PBS solution, and equilibrated in 30% sucrose PBS solution overnight at 4°C. Sagittal sections were cut serially at a 14-μm thickness and then collected on CLIP CRNR Excell slides (Thermo Fisher Scientific, #22-037-247). For immunofluorescence staining, slides were incubated with a primary antibody at 4°C overnight, washed with PBS-Tween, and then incubated with a secondary antibody for 2 hours at room temperature. Antibodies against TH (Pel-Freez, #P40101), HuC (Life Technologies, #A-21271), and activated caspase-3 (BD Biosciences, #559565) were used as primary antibodies. Secondary antibodies were conjugated with Alexa Fluor 488, Alexa Fluor 568, and Alexa Fluor 647 (Life Technologies). 4′,6-Diamidino-2-phenylindole (Life Technologies, #S36973) was used for nuclear staining. Fluorescent images were taken with a Leica SP5X scanning confocal microscope at the Confocal and Light Microscopy core facility at Dana-Farber Cancer Institute.
Zebrafish paraffin sectioning, H&E staining, and immunohistochemistry with primary antibodies against TH (Pel-Freez, #P40101) were performed at the DF/HCC Research Pathology Core.
Cell lines and cell culture
The NGP neuroblastoma cell line was stocked in our laboratory and maintained in RPMI 1640 medium (Gibco) supplemented with 10% fetal bovine serum (Sigma-Aldrich) and 1× penicillin-streptomycin (Gibco). The absence of mycoplasma contamination was confirmed using the MycoAlert Mycoplasma Detection Kit (Lonza, #LT07-418).
CRISPR-Cas9 technology was used to knock out ARID1A in NGP cells. A control gRNA (Ctrl_gRNA, 5′-GAC CGG AAC GAT CTC GCG TA-3′) or each of huARID1A-specific gRNAs (ARID1A_gRNA #1m, 5′-GCT GCT GTG TCG TCG ACT GC-3′ and ARID1A_gRNA #2, 5′-GCC TGC TGG GAG AGC GTC GA-3′; designed by CRISPR Design web tool; http://crispr.mit.edu) was cloned into lentiCRISPR v2 vector (Addgene, #52961) (47, 48). Each of these vectors was cotransfected with pMD2.G and psPAX2 into human embryonic kidney (HEK) 293T cells using FuGENE 6 transfection reagent (Promega, #PRE2693). Cell culture medium containing the lentivirus was filtered by 0.45-μm nitrocellulose membrane and infected with neuroblastoma cells in the presence of polybrene (8 μg/ml). The cells were treated with puromycin (Sigma-Aldrich, #P8833) at 48 hours after infection to eliminate the noninfected cells. Single cell–derived ARID1A homozygous mutant clones were selected by limiting dilution, and the mutations were confirmed by Sanger sequencing.
To generate ARID1A+/− cells, ARID1A_gRNA #2 was cloned into the pSpCas9(BB)-2A-GFP (PX458) vector (Addgene, #48138) (49) and transfected into NGP cells using FuGENE HD transfection reagent (Promega, #PRE2311). Cells were incubated for 72 hours, and GFP-positive single cells were sorted by flow cytometry and genotyped by Sanger sequencing.
Cell growth assay
Cells were seeded in a 96-well plate at a density of 2500 cells per well. Relative cell growth at days 1, 3, 5, and 7 was evaluated by CellTiter-Glo Luminescent Cell Viability Assay (Promega, #G7571).
For cell growth assay with ARID1B knockdown in control, ARID1A+/−, ARID1A−/− #1, and ARID1A−/− #2 NGP cells, a control luciferase shRNA (shLuc, 5′-CGT ACG CGG AAT ACT TCG A-3′) (50) or each of ARID1B-specific shRNAs (shARID1B#1, 5′-GCC GAA TTA CAA ACG CCA TAT-3′ and shARID1B#2, 5′-GGG TTT GGC CCA GGT TAA TAA-3′; https://portals.broadinstitute.org/gpp/public/) was cloned into the CS-RfA-ETV lentiviral expression vector (a gift from H. Miyoshi at RIKEN Center for Computational Science, Japan). Each of the cell line was cotransfected with the doxycycline-inducible shLuc, shARID1B #1, or shARID1B #2, together with pMD2.G and psPAX2 into HEK293T cells using the FuGENE HD transfection reagent (Promega, #PRE2311). Cell culture medium containing the lentivirus was filtered by 0.45-μm nitrocellulose membrane and infected with neuroblastoma cells in the presence of polybrene (8 μg/ml). The infected cells were sorted by the coexpressed EGFP at 72 hours after infection. Cells were pretreated with doxycycline (1.5 μg/ml) for 3 days to induce the knockout of ARID1B before the regular cell growth assay, which is described above.
The half-maximal inhibitory concentration assay
Cells were seeded into a 96-well plate at a density of 10,000 cells per well 1 day before treatment. After 3 days of treatment with vehicle [N,N′-dimethylformamide (DMF)] or cisplatin at various concentrations, relative cell growth was evaluated by CellTiter-Glo Luminescent Cell Viability Assay (Promega, #G7571) as a percentage of the vehicle (DMF) control. The half-maximal inhibitory concentration value was determined with GraphPad Prism 7.01.
Cell transwell invasion and migration assay
A transwell insert with 8.0-μm pore polycarbonate membrane (Corning, #3403) was coated with collagen. Cells were seeded into the upper compartment with serum-free media. Media containing 10% fetal bovine serum were added to the lower compartment. After incubation for 48 hours, the cells were visualized by staining with 0.1% crystal violet. The cells remaining on the upper membrane were removed by cotton wool. The number of migrated cells in five representative fields per membrane was counted. Experiments were performed in triplicate, with statistical significance determined by two-tailed unpaired t test.
Chromatin immunoprecipitation coupled with massively parallel DNA sequencing
ChIP-seq was performed as described previously (35). Human neuroblastoma cells were cross-linked in 11% formaldehyde at room temperature for 15 min, followed by tris-HCl buffer (pH 7.5) quenching. Nuclear extracts were generated with the Nuclei Isolation Kit: Nuclei EZ Prep (Sigma-Aldrich, #NUC101), and chromatin was fragmented by sonication. Anti-BRG1 rabbit monoclonal (EPNCIR111A) antibody (Abcam, #ab110641), anti-SS18 rabbit monoclonal (D6I4Z) antibody (Cell Signaling Technology, #21792), anti-ARID2 mouse monoclonal (E-3) antibody (Santa Cruz Biotechnology, #sc166117X), anti-H3K4me3 rabbit polyclonal antibody (Abcam, #ab8580), anti-H3K27ac rabbit polyclonal antibody (Abcam, #ab4729), and anti-ARID1A mouse monoclonal (PSG3) antibody (EMD Millipore, #04-080) coupled with Dynabeads M-270 Epoxy (Life Technologies, #14301) were used for immunoprecipitation. ChIP DNA was purified and subjected to library preparation for Illumina NextSeq 500 Next-Generation Sequencing at the Dana-Farber Cancer institute Molecular Biology Core Facilities. ChIP-seq data are available from Gene Expression Omnibus (GEO) (series GSE134632).
ChIP-seq data analysis and display
ChIP-seq reads were aligned to human reference genome (hg19) using Bowtie2/2.2.9 (51) with default parameters. SAMtools/1.3.1 (52) was used to remove duplicated reads and reads with low mapping quality (mapQ < 20). SS18 and ARID2 peaks were called against input reads using MACS2/2.1.1 (53) with default parameters. BRG1 and H3K4me3 peaks were called against input reads using MACS2/2.1.1 at q = 1 × 10−3. H3K27ac peaks were called twice using MACS (53) with parameter sets –keep-dup=auto –p 1e-9 and –keep-dup=all –p 1e-9 and then collapsed together for downstream analysis. Peaks overlapping with MYCN amplified regions (chr2: 16004820-17229328) and ENCODE blacklisted regions were discarded for downstream analysis.
An SS18 peak with at least a 50% overlap with a BRG1 peak was identified as a BAF binding site, while an ARID2 peak with at least a 50% overlap with a BRG1 peak was identified as a PBAF binding site. BAF or PBAF binding sites were determined in each condition (control, ARID1A−/− #1 and ARID1A−/− #2). The BAF and PBAF binding sites in control NGP cells were used for annotation with R Bioconductor package ChIPseeker (31). For this annotation, genes with the closest TSS were considered to be BAF/PBAF complex target genes.
The regions of BAF or PBAF binding sites in all three conditions were then merged together by bedtools 2.26.0 (54) as BAF or PBAF binding sites. To define the alterations of BAF/PBAF binding upon ARID1A loss, ChIP-seq read counts were generated by bedtools at each site, and read counts for each binding site were normalized to per million mapped reads (RPM). Input RPM values for each BAF/PBAF binding site in each condition were subtracted from SS18/ARID2 RPM values in that condition, and all the negative values gained at this step were set to 0. Fold changes were evaluated at each BAF/PBAF binding site using the normalized RPM values mentioned above, with a pseudocount of 0.1. If the fold change of BAF/PBAF binding at a site was greater than 1.5-fold in both ARID1A mutants, then the site was called a strengthened BAF/PBAF binding site. If, on the other hand, the fold change was less than two/threefold in both ARID1A mutants, the site was called a weakened BAF/PBAF binding site. Other sites were designated stable BAF/PBAF binding sites.
TSS-proximal/distal sites were classified as previously described (16). Active TSS was defined as all the TSS annotated in the UCSC hg19 gene track that overlap with H3K4me3 peaks in any of the three conditions. BAF/PBAF binding sites overlapping with both an H3K4me3 peak and an active TSS were designated TSS proximal. Those more than 2 kb away from an active TSS and more than 1 kb away from an H3K4me3 peak were considered TSS distal. Other ambiguous sites were removed for TSS-proximal/distal studies. H3K27ac peaks in all the conditions were merged together. Similar to BAF/PBAF binding sites, H3K27ac peaks overlapping with both an H3K4me3 peak and an active TSS were designated promoters, while those more than 2 kb away from an active TSS and more than 1 kb away from an H3K4me3 peak were considered enhancers.
Bigwig files for display were made RPM using MACS2/2.1.1 with input subtracted. Heat maps for BAF/PBAF binding were generated by deepTools 3.0.2 (55). Every line of the heat maps shows the input-subtracted values with negative values set to 0.
RNA-seq data analysis
RNA-seq reads were aligned to hg19 using HISAT2 2.1.0 (56) with default parameters. GFold v1.1.2 (57) was used to generate reads per kilobase of exon per million reads mapped (RPKM) values and calculate log2 fold change between ARID1A mutants and control cells from RNA-seq data. RNA-seq data are available from GEO (series GSE134632).
Each TSS-distal BAF or PBAF binding site was annotated to the gene with the closest active TSS and classified as strengthened, stable, or weakened site as defined above by ChIP-seq signals. The log2 fold change between ARID1A mutants and control cells was plotted for each category. Statistical significance was determined with two-sided Mann-Whitney test (http://astatsa.com/WilcoxonTest/).
GSEA and GO analysis
Normalized expression values for individual samples of RNA-seq were obtained from DESeq2 (58) using the variance-stabilizing transformation of the raw counts and were used for GSEA. Genes from whole transcriptomic data were ranked on the basis of the stat from DESeq2 in ARID1A mutant versus control NGP cells. The preranked option of GSEA (59) was run with 2500 permutations for statistical evaluation. GSEA was performed with signatures from version 6.0 of the Molecular Signatures Database (www.broadinstitute.org/gsea/msigdb/index.jsp). A generic gene signature for mesenchymal tumor (38) was also used as a gene set for GSEA. The GO analysis was performed using the Enrichr program (http://amp.pharm.mssm.edu/Enrichr/) (32), and the top 10 GO terms ranked by combined score were shown.
The genomic sequence within 100 base pairs of the center of each weakened BAF binding site in ARID1A mutant NGP cells was generated from the UCSC website with repeats masked with “N.” Motif analysis was performed by MEME-ChIP v4.12.0 (http://meme-suite.org/tools/meme-chip) (60).
To analyze the hyperplastic, apoptotic, and proliferative status of the sympathoadrenal cells before tumor onset, we scanned all sections from each fish by confocal z-stack. A single representative section containing the largest number of sympathoadrenal cells or activated caspase-3–positive sympathoadrenal cells or the largest percentage of EdU+ sympathoadrenal cells in the IRG was selected and quantified for each fish. Data were compared with either the two-tailed unpaired Welch t test (hyperplasia analysis) or the two-tailed Fisher’s exact test (apoptosis analysis).
To analyze the proliferative capacity of the sympathoadrenal cells before tumor onset, all sections from each individual fish were scanned by confocal z-stack. A single representative section containing the largest percentage of EdU+ sympathoadrenal cells to the total number of sympathoadrenal cells in the IRG was selected and quantified for each individual fish. The data were analyzed by two-tailed unpaired Welch t test.
Acknowledgments: We thank J. R. Gilbert for review of the manuscript and comments; G. Thurston, D. Debiasi, M. T. Coute, and M. Alves for care of the zebrafish; Z. Herbert of the Dana-Farber Molecular Biology Core Facility for genomics support; and Y. Zhou and S. Yang of Boston Children’s Hospital for comments. Funding: This work was supported by an NIH grant R01-CA180692 (A.T.L. and J. M. Maris co-PI grant). T.T. was supported by a research grant from Children’s Neuroblastoma Cancer Foundation, a Friends of DFCI award, and a Friends for Life Fellowship from Dana-Farber Cancer Institute. H.S. was supported by a Young Investigator Grant from Alex’s Lemonade Stand Foundation, a Rally Foundation for Childhood Cancer Research, and a Friends for Life Neuroblastoma Fellowship from Dana-Farber Cancer Institute. M.W.Z. was supported by a Young Investigator Grant from Alex’s Lemonade Stand Foundation, a postdoctoral fellowship from the Charles A. King Trust, and the Claudia Adams Barr Innovative Basic Science research program. A.D.D. was supported by funding from the Damon-Runyon Cancer Research Foundation (DRSG-24-18), the Alex’s Lemonade Stand Foundation, and the CureSearch for Children’s Cancer Foundation. Author contributions: H.S., T.T., C.K., and A.T.L. designed the study and experiments. H.S. and T.T. performed the experiments. H.S., T.T., and B.J.A. performed the bioinformatics analysis. H.S., T.T., B.J.A., A.D.D., M.W.Z., C.K., and A.T.L. analyzed the data. H.S., T.T., and A.T.L. drafted and finalized the manuscript with inputs from all other authors. Competing interests: C.K. is the scientific founder, a fiduciary board of directors member, a scientific advisory board member, a shareholder, and a consultant for Foghorn Therapeutics Inc. (Cambridge, MA). B.J.A is a shareholder in Syros Pharmaceuticals Inc. The other authors declare that they have no potential conflicts of interest. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. ChIP-seq and RNA-seq data are available from GEO (series GSE134632). Additional data related to this paper may be requested from the authors.