Chromatin remodeling protein BPTF mediates chromatin accessibility at gene promoters in planarian stem cells

Abstract

Background

The regulation of chromatin accessibility is essential in eukaryotic cells as one of several mechanisms that ensure gene activation occurs at appropriate times and in appropriate cell types. Accordingly, mutations in chromatin remodeling proteins are linked to many different developmental disorders and cancers. One example of a chromatin protein that has been linked to both developmental abnormalities and cancer is BPTF/NURF301, the largest subunit of the Nucleosome Remodeling Factor (NuRF) complex. The BPTF subunit is not only important for the formation of NuRF but also helps direct its activity to particular regions of chromatin by preferentially binding histone H3 lysine four trimethylation (H3K4me3). Notably, defects caused by knockdown of bptf in Xenopus embryos mimic those caused by knockdown of wdr5, a core subunit of all H3K4me3 methyltransferase complexes. However, the mechanistic details of how and where BPTF/NuRF is recruited to regulate gene expression vary between studies and have been largely tested in vitro and/or in cultured cells. Improving our understanding of how this chromatin remodeling complex targets specific gene loci and regulates their expression in an organismal context will provide important insight into how pathogenic mutations disrupt its normal, in vivo, cellular functions.

Results

Here, we report our findings on the role of BPTF in maintaining chromatin accessibility and essential function in planarian (Schmidtea mediterranea) stem cells. We find that depletion of planarian BPTF primarily affects accessibility at gene promoters near transcription start sites (TSSs). BPTF-dependent loss of accessibility did not correlate with decreased gene expression when we considered all affected loci. However, we found that genes marked by Set1-dependent H3K4me3, but not MLL1/2-dependent H3K4me3, showed increased sensitivity to the loss of BPTF-dependent accessibility. In addition, knockdown of bptf (Smed-bptf) produces loss-of-function phenotypes similar to those caused by knockdown of Smed-set1.

Conclusions

The S.mediterranea homolog of NuRF protein BPTF (SMED-BPTF) is essential for normal homeostasis in planarian tissues, potentially through its role in maintaining chromatin accessibility at a specific subset of gene promoters in planarian stem cells. By identifying loci that lose both chromatin accessibility and gene expression after depletion of BPTF, we have identified a cohort of genes that may have important functions in stem cell biology.

Background

Planarians are a robust and powerful model for studying in vivo stem cell function. These animals maintain a large population of adult stem cells that they use for both unlimited tissue turnover and to regenerate new organs and tissues after major injury [1]. Unlike most cultured stem cells grown in optimized conditions, the planarian stem cell population is a highly heterogeneous mix of multi and pluripotent cells (also known as neoblasts) that proliferate and differentiate asynchronously in response to injury signals [2]. Unsurprisingly, the planarian stem cell population is also transcriptionally heterogeneous [3,4,5,6,7,8,9] with some stem cells expressing various lineage-specific genes (e.g., myoD or foxA) in addition to common stem cell markers (e.g., piwi-1). Studies have shown that these lineage-specific stem cells are required for the maintenance and regeneration of their corresponding tissues [10,11,12]. Additional research has revealed that positional cues [13], cell cycle phase [14], mitochondrial state [15], and specific signaling pathways [16] are all important regulators of stem cell heterogeneity and function. However, the details about how such changes in cell state and cell signaling are translated into transcriptional changes remain unclear.

One category of proteins with the potential to act as major regulators of transcriptional dynamics in planarian stem cells is chromatin modifying proteins. This group of proteins includes enzymes that add covalent modifications to histone proteins, those that remove them, and complexes that regulate (or “remodel”) chromatin accessibility by shifting and/or removing nucleosomes [17,18,19]. These modifiers each play important roles in the multilayered, synergistic mechanisms that regulate gene expression in eukaryotic cells. Several studies have identified the planarian homologs of various chromatin proteins and uncovered loss-of-function phenotypes for many of them [20,21,22,23,24,25,26,27,28]. More recently, as new protocols have been developed, better genomes have been assembled, and the sensitivity of chromatin assays has increased, researchers in the planarian field have made substantial progress in understanding the chromatin state and how it is regulated in planarian cells [25, 26, 28,29,30,31,32]. Together, these studies have revealed both predictable and unexpected findings about the chromatin states of planarian cells and how they are regulated, suggesting that further research in this model will uncover important mechanistic details about how chromatin mechanisms regulate both regeneration and stem cell function.

In this study, we examined both the functional and molecular roles of the chromatin remodeling protein BPTF in planarian stem cells. BPTF is the largest subunit of the Nucleosome Remodeling Factor (NuRF), an ISWI-containing chromatin remodeling complex that executes its function by shifting nucleosomes bidirectionally in 10 bp increments [33,34,35]. BPTF is not only an essential subunit of NuRF but also plays a major role in focusing its activity at specific genomic loci through its PHD2 domain, which preferentially binds H3K4me3 [36, 37]. BPTF binding is not only highly specific for H3K4me3 (over other methylated lysines in the H3 tail or other states of H3K4 methylation), but knockdown of bptf in Xenopus embryos also phenocopies the axial, blood, and gut defects seen in embryos with reduced H3K4me3 (via wdr5 knockdown) [37]. Together, these data support the hypothesis that BPTF mediates chromatin remodeling and downstream functional effects at genomic loci with specific chromatin features.

In keeping with previous studies that uncovered essential roles for BPTF in mammalian stem cells [38,39,40], we found that knockdown of Smed-bptf by RNA interference (RNAi) caused defects in planarian stem cell function and tissue homeostasis. To identify the specific gene loci affected by the loss of BPTF, we used the Assay for Transposase-Accessible Chromatin with sequencing (ATAC-seq), Chromatin Immunoprecipitation plus sequencing (ChIP-seq), and RNA sequencing (RNA-seq) to identify changes in both the chromatin state and transcriptional output in stem cells isolated from bptf(RNAi) planarians. These data revealed that BPTF is important for maintaining chromatin accessibility at the promoters of over 3000 genes, although only a subset of these loci showed corresponding changes in gene expression. We found that genes with both H3K4me3 and ATAC-seq peaks at their promoters had robust and consistent patterns of chromatin accessibility focused near Transcription Start Sites (TSSs), whereas genes with ATAC-seq peaks only (no H3K4me3) were highly accessible but with less consistent patterns of open chromatin. These data support the model in which H3K4me3 recruits NuRF activity to specific loci via BPTF binding.

The S.mediterranea genome encodes two H3K4me3 KMTases, Set1 and MLL1/2, and RNAi of these enzymes causes distinct phenotypes: knockdown of Smed-set1 causes progressive stem cell dysfunction, tissue regression, and death, whereas knockdown of Smed-mll1/2 induces motility defects due to the gradual loss of epidermal cilia [23, 29]. Importantly, these phenotypes are clearly linked to their distinct genomic targets in planarian stem cells i.e., genes associated with stem cell function and ciliogenesis respectively [29]. Notably, the changes we detected in chromatin accessibility in bptf(RNAi) stem cells were more likely to correlate with changes in gene expression if their promoters were marked by Set1-dependent H3K4me3. This finding aligns with the strong similarity in set1(RNAi) and bptf(RNAi), but not mll1/2(RNAi), phenotypes. These data strongly support conservation in function across organisms, as loss of H3K4me3 and loss of BPTF also produce highly similar phenotypes in Xenopus. Additionally, our findings identify those H3K4me3-marked genes at which BPTF binding is critical for their stable expression. Given that BPTF has been shown to be essential for stem cell function in both mammals and planarians, our results may point to specific genes and mechanisms that are fundamental to stem cell biology.

Methods

Homolog identification and phylogenetic tree construction

To identify putative BPTF orthologs in the planarian species Schmidtea mediterranea, we used validated human (Q12830), mouse (A2A654), and Drosophila melanogaster (Q9W0T1.2) BPTF/NURF301 protein sequences in the UniProtKB/Swiss-Prot database [41] as query sequences for protein BLAST (BLASTP) searches against the S2F19 HAP1 PRJNA885486 genome assembly [31] using an evalue cutoff 1.0e-10 (accessed at WormBase ParaSite [42]). The top hit for all queries was h1SMcG0002696. Reciprocal Best Hits (RBH) [43] was used to validate this putative ortholog against all query sequences. Other BPTF orthologs for phylogenetic tree construction were identified from the ENSMBL database [44] and the associated protein sequence downloaded from UniProtKB/Swiss-Prot [41]. Species chosen for inclusion were model organisms (Danio rerio), comparator species (Crassostrea gigas, Ixodes scapularis, Ciona intestinalis) and other invertebrates based on their phylogenetic position (Hofstenia miamia and Macrostomum lignano). The sequence of the conserved PHD2 domain was identified and isolated from each full-length predicted protein using SMART sequence analysis tool [45]. Both full-length proteins and the isolated PHD2 domains were each aligned by Mafft [46] using default settings and these alignments were subsequently curated using the BMGE software [47]. Phylogenetic trees were constructed in the https://ngphylogeny.fr/ interface [48] using PhyML with default settings except best of NNI and SPR was used for tree topology search and approximate Bayes was used for branch support [49]. Trees were visualized and re-rooted using the iTOL online tool [50]. Human p300 (Q09472) was used as an outgroup to root both the full-length and isolated PHD-domain trees.

Gene cloning and RNAi

Two separate, non-overlapping regions of the identified Smed-bptf gene (h1SMcG0002696) were cloned from a cDNA library using the following gene-specific primers: bptf-1_F = TGAGTAATGTCAATATGAAACC, bptf-1_R = AATAGTCCACATCCGTATATCT, bptf-2_F = AATGTGATAAAACACACGATAC, bptf-2_R = GAAGTTCAGTTAAAAAGTAGGC. We then used established methods [51] to clone each of these sequences into a T4P vector with flanking T7 binding motifs, transform the resulting plasmids into the bacterial strain Ht115 [52], induce these bacteria to express the cloned sequence as dsRNA (1mM final IPTG, 2 h shaking at 37 C), pellet, rinse, and mix the resulting bacterial pellet with calf liver paste. We then fed animals with this paste as indicated in the results, using a previously generated T4P vector containing a segment of the Caenorhabditis elegans unc-22 gene as a control [53]. Each non-overlapping bptf-RNAi construct (1 & 2) induced the same phenotype in multiple independent experiments, after which the bptf-1 RNAi construct was used for all genomic/transcriptomic experiments.

Planarian culture and radiation

An asexual strain (CIW4) of S.mediterranea worms were maintained as described [54] with gentamicin supplemented at 50ug/ml. Worms used in experiments were starved ≥ 7 days prior to selection for experiments. Irradiation was done with a Gammacell 40 Exactor (MDS Nordion) at a dose rate of 85 Rad/minute.

Live worm imaging and whole-mount in situ hybridization (WISH)

Live worms images were captured using a Leica M205-FA. Riboprobe for piwi-1 [55] was made by in vitro transcription with T7 polymerase from the pPR-T4P vector as previously described using DIG labeling [56]. Whole-mount in situ hybridization (WISH), including tyramide amplification of the fluorescent piwi-1 riboprobe, was performed as previously described [57]. Worms were then mounted in modified ScaleA2 (30% glycerol/4 M urea/0.1% TritonX-100/2.5% DABCO) and imaged using a confocal microscope (LSM 700, Zeiss). Quantitation of piwi-1⁺ cells in the confocal images was performed using ImageJ [58]. Briefly, we used the “analyze particles” function to count particles (i.e., piwi-1⁺ cells) that met defined size and circularity criteria (10–120 pixels, 0–1 circularity). We then averaged the number of piwi-1⁺ cells across the animals in each RNAi condition (n≥9).

Flow cytometry

Animals were first minced with a clean razor blade, then dissociated in cold calcium- and magnesium-free buffer with 1% BSA (CMFB) over the course of 15 min with occasional pipetting using a P1000 filter top. Dissociated worms were then filtered through 100 μm followed by 40 μm cell strainers to remove non-dissociated tissue. Remaining cells were pelleted by low-speed centrifugation at 300rcf (5 min) and resuspended in room temperature CMFB + Hoechst 33,342 (10 µg/ml). Dissociated cells were incubated with Hoechst 33,342 at room temperature in the dark for 1.5 h. Propidium iodide was added (5 µg/ml) for the last 5 min to assess cell viability. Viable > 2n “X1” stem cells [59] were then collected using a MoFlo Legacy (Beckman Coulter) or BD Influx (BD Biosciences). Both control (unc-22) and bptf(RNAi) animals were processed in parallel, using both available FACS machines in the Stowers Institute Core facility concurrently to minimize variation between control and bptf(RNAi) cells. Cells were stored on ice briefly until all replicates were isolated and then processed immediately (see below). Post-sort analysis was conducted with X1-gated cells to ensure purity before proceeding with downstream applications.

Assay for transposase-accessible chromatin with sequencing (ATAC-seq)

Two replicates of 50 K stem cells per control(RNAi) and bptf(RNAi) worms were isolated by flow cytometry (see above) and subjected to ATAC according to the original published protocol [60] and using Tn5 transposase enzyme purchased from Illumina as part of the Nextera DNA library prep kit (#FC-121-1030). A total of ten amplification cycles were used in the generation of these ATAC-seq libraries, a number determined by the mid-amplification qPCR assessment recommended in the protocol [60]. The resulting ATAC-seq libraries were sequenced on the Illumina NextSeq 500 (paired-end reads) and quality check (QC) was done using fastqc (version 0.11.9) [61]. Fragment size distributions were generated using ATACseqQC (version 1.26.0) with the fragSizeDist function [62]. Bowtie2 (version 2.4.4) [63] was used to align sample reads to the S3h1 S.mediterranea genome [31] using the parameters --very-sensitive -X 1000. After alignment, samtools [64] was used to create BAM files for each sample with those alignments that are properly paired (samtools view -h -b -f 3). Samples were then subsampled based on library complexity and processed as previously described [65]. PCR duplicated reads were removed using the MarkDuplicates function in Picard (version 2.26.4) (http://broadinstitute.github.io/picard).

To compensate for Tn5 insertion, a Tn5 shift (+ 4 and − 5 bp) was applied using the bedpeTn5shift.sh script [65]. ATAC-seq peaks were then called using the MACS2 [66] callpeak function (version 2.2.7.1) and parameters -g 8.40e8 -q 0.05 --keep-dup all. Consensus peaks present in both replicates were identified using BEDtools [67] and the annotation of MACS2 peaks to gene models was done using ChIPSeeker (version 1.38.0) [68, 69]. deepTools (version 3.5.1) [70] was used to generate normalized coverage files (bigWig files) and profile plots, using parameters bamCoverage --binSize 5 --smoothLength 60 --normalizeUsing RPKM. The Integrative Genomics Viewer (IGV version 2.16.2) was used to visualize bigWig files [71]. The identification of differentially accessible ATAC-seq peaks was done using csaw (version 1.36.0) [72] as recommended in a previous study [65]. Further quantitative and statistical analysis was performed using R (version 4.3.1) and figures were created using ggplot2 [73].

Gene expression analysis

Four replicates of 100 K X1 stem cells per control(RNAi) and bptf(RNAi) condition per time point were isolated by flow cytometry (see above) and collected directly into Trizol LS regent (a version of Trizol optimized for liquid samples, cat #10296010). The isolated cells were then homogenized and total RNA extracted per the Trizol manufacturer protocol. RNAseq libraries were then generated using the Illumina TruSeq kit and sequenced in 50 bp single reads using the Illumina HiSeq 2500. Single-end reads were then aligned to the S3h1 S.mediterranea genome using Hisat2 (version 2.2.1) [31, 74] and transcript quantification was performed using StringTie (version 2.1.7). Differential gene expression analysis was then performed using DESeq2 (version 1.42.0) [75, 76]. Genes were considered differentially expressed if they had p-adj < 0.05. Further quantitative and statistical analysis was performed using R (version 4.3.1) and figures were created using ggplot2 [73]. For correlation analyses, regression lines were fitted using the linear model (lm method) in geom_smooth() from ggplot2.

Chromatin-immunoprecipitation plus DNA sequencing (ChIP-seq)

Chromatin-immunoprecipitation was performed as previously described [29, 77]. For each replicate/condition, 300 K planarian X1 stem cells were isolated by FACS (see above), cross-linked in 4% PFA, and mixed with 10 M (cross-linked) Drosophila S2 cells before proceeding to the next steps. Chromatin was sheared in the S220 Focused Ultrasonicator (Covaris) with silia beads and immunoprecipitation was performed with an antibody to H3K4me3 (Cell Signaling #9751). Isolated genomic DNA from both input and ChIP samples were used to generate libraries using the Illumina TruSeq kit and then sequenced (single end) on the Illumina HiSeq 2500. Quality checks of all resulting reads (QC) were done using fastqc (version 0.11.9) [61] followed by alignment to the S3h1 S.mediterranea genome with bowtie2 (version 2.4.4) [76, 78]. PCR duplicated reads were removed using MarkDuplicates function in Picard (version 2.26.4) (http://broadinstitute.github.io/picard). Peak calling was done using the MACS2 [66] callpeak function (version 2.2.7.1) with–g 8.40e8 --nomodel and -q 0.01 parameters. Differential analysis between set1(RNAi), mll1/2(RNAi), and bptf(RNAi) H3K4me3 ChIP libraries with their individual matching control(RNAi) libraries was performed using diffReps (version 1.55.6), paramaters --btr --bco --nsd --frag 250 [79]. Consensus peaks present in both replicates were identified using BEDtools [67] and peaks annotation to gene models was done using ChIPSeeker (version 1.38.0) [68, 69]. deepTools (version 3.5.1) [70] was used to generate normalized coverage files (bigWig files) and for creating profile plots (bamCoverage --normalizeUsing RPKM). Integrative Genomics Viewer (IGV version 2.16.2) was used to visualize bigWig files [71]. The GO (gene ontology) enrichment analysis was done using enricher function from package clusterProfiler (version 4.10.0) [80, 81]. Further quantitative and statistical analysis was performed using R (version 4.3.1) and figures were created using ggplot2 [73].

Results

The genome of Schmidtea mediterranea encodes a homolog of NuRF protein BPTF

BPTF is a H3K4me3-binding protein and the largest subunit of the NuRF chromatin remodeling complex (Fig. 1A). To characterize its functional and molecular role in planarian stem cells, we first identified the likely homolog of bptf in the S.mediterranea genome using a reciprocal BLAST strategy [43]. This analysis identified a single gene encoding a protein with 27% sequence identity to human BPTF and 28% identity to the Drosophila homolog NURF301 (calculated using Clustal Omega; DmNURF301 has 34.5% identity with human BPTF in the same matrix). Despite the modest conservation of overall sequence between these three BPTF homologs, their predicted protein domains are highly conserved (Fig. 1B). Moreover, after aligning their PHD2 motifs and examining these sequences more closely, we found that all four aromatic residues essential for H3K4me3 binding [82] are conserved across human, Drosophila,

Fig. 1

The genome of Schmidtea mediterranea encodes a homolog of BPTF, the H3K4me3-binding subunit of the Nucleosome Remodeling Factor (NuRF) complex. A) Cartoon showing the NuRF (Nucleosome Remodeling Factor complex), a member of the ISWI complex family that uniquely contains the H3K4me3-binding subunit BPTF. Orange circles represent methyl groups (in set of three = trimethylation), asterisks represent the key amino acids required for H3K4me3 binding. B) Illustration comparing the domain structures of human BPTF, Drosophila NURF301 and a BPTF/NURF homolog identified in the planarian species Schmidtea mediterranea (Smed). C) Alignment of the PHD2 domain from human (NP_872579.2), Drosophila, and SMED-BPTF. Those amino acids known to be critical for its H3K4me3-binding [36] are highlighted and asterisked. D) Phylogenic tree constructed using complete BPTF protein sequences for ten species. Numbers indicate Bayesian posterior probabilities. E) Phylogenic tree constructed using the PHD2 domain of each of the ten BPTF protein sequences in D. Numbers indicate Bayesian posterior probabilities

and SMED-BPTF (Fig. 1C, asterisks). To address the question of SMED-BPTF homology further, we used this sequence and the BPTF homologs of nine other species to construct phylogenetic trees (Fig. 1D, E). Trees constructed with either full-length BPTF sequences (Fig. 1D) or their PHD2 motifs only (Fig. 1E) suggested that human BPTF, Drosophila NURF301 and SMED-BPTF are orthologous to each other.

Knockdown of chromatin remodeling protein BPTF mimics knockdown of KMTase Set1

Having identified a homolog of BPTF in S.mediterranea, we then asked if it is required for in vivo function. To address this question, we used RNAi to knock down Smed-bptf in adult planarian worms. In parallel, we also knocked down Smed-set1, Smed-mll1/2, and a control gene (C. elegans unc-22) using previously reported RNAi constructs [29] for comparison. As shown in previous studies [23, 29], set1(RNAi) animals showed defects in tissue homeostasis, including progressive head regression (Fig. 2A), within two weeks of dsRNA treatment and were all dead within three weeks (Fig. 2B). In contrast, mll1/2(RNAi) worms did not show signs of head regression or death even after several weeks (Fig. 2A, B) but did develop motility defects due to loss of ventral cilia [23, 83]. RNAi of Smed-bptf caused a penetrant and morphologically similar phenotype to that of set1(RNAi), although its progression developed more slowly (Fig. 2A, B). Notably, this was true for two separate, non-overlapping bptf RNAi constructs (Fig. 2B), supporting the conclusion that this phenotype is specific to Smed-bptf knockdown and not an off-target effect. However, we note that because the tissue regression phenotypes of set1(RNAi) and bptf(RNAi) are relatively quick (2–3 weeks) compared to the slower onset of the mll1/2(RNAi) motility defect (3–4 weeks), we cannot conclude from these data that SMED-BPTF preferentially operates at Set1 target genes.

Fig. 2

RNAi knockdown of Smed-bptf phenocopies that of Smed-set1, but not Smed-mll1/2. A) Live images showing the morphological phenotypes of bptf(RNAi) planarians (S.mediterranea) compared to control(RNAi), set1(RNAi), and mll1/2(RNAi) animals. Scale bars = 200 μm. B) Survival plot of set1(RNAi) and bptf(RNAi) animals compared to control RNAi worms. Two non-overlapping constructs are shown for each gene, set1 and bptf, to control for off target effects. C) Whole mount In situ hybridization (WISH) in wild-type worm using a riboprobe to bptf. D) Plots generated from previously published single cell RNA-seq data [7] using the publicly available Shiny App https://simrcompbio.shinyapps.io/bbp_app/ to show bptf expression in specific cell types. E) Schematic of the experimental setup used to test planarian stem cell function in RNAi worms. In normal conditions, a small but significant number of planarian stem cells will survive 1250 rad (12.5 Gy) γ-radiation (3 days post-irradiation, or dpi), then resume proliferating (7 dpi) to restore the population [24, 26, 84]. F) Survival plot of set1(RNAi), mll1/2(RNAi), bptf(RNAi) and control(RNAi) animals with (dotted line) and without (solid line) radiation treatment. G) Representative Fluorescence In situ hybridization (FISH) images of control(RNAi) and bptf(RNAi) planarians stained with the stem cell marker piwi-1. Scale bars = 250 μm. H) Quantitation of piwi-1+ cells per animal for all animals included in E. For each condition (RNAi treatment, +/- radiation, time point) n = 9–12. Statistical significance determined using student’s t-test, *** = p-value < 0.001

Head and tissue regression in planarians can be caused by the failure of many different cellular processes, including those that are critical to planarian stem cells (as these are the only cells that divide and replenish planarian tissues). To assess whether the expression pattern of Smed-bptf showed enrichment in a particular tissue or cell type, we performed in situ hybridization using a riboprobe for Smed-bptf (Fig. 2C). Smed-bptf was detected in a broad pattern of expression with some darker staining outlining the cephalic ganglia and the stem cell dense “tail stripe” (Fig. 2C). Published RNA-sequencing datasets generated from FACS-isolated stem cells [85] and single cell RNA-sequencing (scRNA-seq; Fig. 2D) [7] both indicate that Smed-bptf expression is moderately enriched in planarian stem cells. However, it is still possible that the bptf(RNAi) phenotype is due to its function in other cells or a combination of cell types.

Previous studies from our lab and others [23, 29, 86] have shown that set1(RNAi) worms show defects in planarian stem cell function, including an inability to recover from low doses of radiation. To assess more specifically if knockdown of Smed-bptf impacts stem cells, we performed a previously established assay of stem cell function [24, 26, 84](Fig. 2E): first, we depleted specific gene transcripts (bptf, set1, mll1/2, and a control gene) with RNAi. Next, we divided each group of RNAi worms into two groups: one that we treated with 1250 rads (12.5 Gy) ionizing radiation and a second that was not irradiated (0 rad). As shown previously [24, 26, 84], this dose of radiation reduces the number of stem cells in each animal, challenging those remaining to repopulate and maintain tissue homeostasis. On days 3 and 7 post radiation, we fixed a subset of bptf(RNAi) and control(RNAi) animals and then monitored the remaining (live) worms for several weeks to score survival and phenotype (Fig. 2E, F). Notably, both set1(RNAi) and bptf(RNAi) animals, but not mll1/2(RNAi) or control animals, displayed an accelerated progression of tissue regression and death after 1250 rad treatment compared to their unirradiated (0 rad) counterparts (Fig. 2F).

We then stained those animals fixed on days 3 and 7 post-radiation with an in situ marker of planarian stem cells (piwi-1; Fig. 2G) and compared the number of stem cells in each subgroup i.e., control(RNAi) +/- radiation, bptf(RNAi) +/- radiation. As expected based on prior experiments [24, 26, 84], control(RNAi) worms treated with 1250 rads showed an initial depletion of most piwi-1 + stem cells, followed by significant restoration of the population by day 7 (Fig. 2G, H). In contrast, bptf(RNAi) worms had significantly fewer stem cells compared to control at 3 days post-radiation treatment and failed to show significant recovery at day 7 (Fig. 2G, H). Notably, set1(RNAi) stem cells also show the same loss-of-function effects in this radiation assay [29], suggesting that the set1(RNAi) and bptf(RNAi) phenotypes may be due to disruption of similar cellular and/or molecular mechanisms.

Knockdown of Bptf leads to loss of chromatin accessibility at gene promoters

After observing that knockdown of Smed-bptf caused significant stem cell dysfunction, we then asked if bptf(RNAi) causes loss of chromatin accessibility at specific loci in these stem cells due to its predicted association with NuRF. We tested this prediction by performing ATAC-seq (60) on isolated X1 stem cells (Fig. 3A) and aligning the data to the S3h1 S.mediterranea genome assembly [31] (Supplemental Figure S1A). We detected thousands of ATAC-seq peaks in both samples, including those with measurable differences between control(RNAi) and bptf(RNAi) and those with no significant change (Fig. 3B). After mapping these peaks to the nearest annotated gene model [68, 69], we found that most ATAC-seq peaks in control(RNAi) stem cells map near promoters and a significant number map to distal intergenic regions (Fig. 3C, Supplemental Figure S1B). We then used csaw [87] to identify those peaks with significant changes in ATAC-seq signal in bptf(RNAi) stem cells.

Fig. 3

Knockdown of bptf leads to loss of chromatin accessibility at gene promoters in planarian stem cells. A) Schematic of the experimental setup used to assay the chromatin state of planarian stem cells isolated from control(RNAi) and bptf(RNAi) animals. B) A representative locus at which there are multiple ATAC-seq peaks, including one that shows significant loss of accessibility after bptf(RNAi). Gene nkx2.2 = h1SMcG0021724. C) Pie chart summarizing the locations of ATAC-seq peaks across the planarian genome (relative to their nearest gene models) in control(RNAi) stem cells. D) Pie chart summarizing the locations of differentially accessible (DA) peaks in bptf(RNAi) stem cells (compared to control(RNAi) stem cells). E) Histogram summarizing the log2FC of all differentially accessible (DA) ATAC-seq peaks in bptf(RNAi) stem cells (binwidth = 0.1 log2FC). 4652 peaks have decreased accessibility (blue bars), 1411 have increased accessibility (orange bars). F) Profile plot comparing ATAC-seq signal in control(RNAi) and bptf(RNAi) stem cells averaged across all gene models and centered at transcription start site (TSS). Signal is represented as Reads Per Kilobase per Million mapped reads (RPKM). The shaded area represents the standard error. G) Profile plot comparing ATAC-seq signal at distal intergenic peaks (identified in C). The shaded area represents the standard error

More than half of the differentially-accessible (DA) peaks we detected in bptf(RNAi) stem cells mapped near gene promoters, with a relatively smaller fraction mapping to distal intergenic regions (Fig. 3D; Supplemental Figure S1C). The majority of these DA peaks lose ATAC-seq signal versus gain (Fig. 3E, Supplemental Figure S1C), as predicted if SMED-BPTF functions in chromatin remodeling via NuRF. This difference was also unambiguously observed when we averaged ATAC-seq signal across all annotated gene models and centered it at Transcription Start Sites (TSSs): chromatin accessibility was significantly lower in bptf(RNAi) stem cells +/- 1 kb of the TSS (Fig. 3F, Supplemental Figure S2A). This was also the case when we compared ATAC-seq signal at defined promoter peaks (Supplemental S2E, F). However, chromatin accessibility was not detectably changed at distal intergenic loci (Fig. 3G, Supplemental Figure S2B) or peaks at “other introns” (Supplemental Figure S2G, H) in bptf(RNAi) stem cells. Additionally, we did not find significant overlap between these BPTF-dependent DA peaks and those previously shown to be differentially accessible after knockdown of the BAF chromatin remodeling complex (also in X1 stem cells [28]). Together, these data support the hypothesis that SMED-BPTF is a conserved subunit of the NuRF complex that helps direct NuRF activity to specific genomic loci in planarian stem cells. However, further biochemical studies will be needed to confirm this hypothesis more conclusively.

Knockdown of Bptf leads to significant changes in gene expression

As changes in chromatin accessibility have the potential to affect transcription, we then asked if knockdown of bptf also impacted gene expression in planarian stem cells. We again isolated X1 stem cells by flow cytometry from control(RNAi) and bptf(RNAi) animals and then extracted their total RNA for RNA-sequencing. We isolated stem cells at two time points post-RNAi treatment (day 9 and day 12 post first RNAi feeding) to assess the timing and consistency of gene expression changes after loss of chromatin accessibility. Indeed, we detected many significant changes at both time points (Fig. 4A), although overall the changes in gene expression were consistent between time points for both up- and down-regulated genes (Fig. 4B, Supplemental Figure S3A). Interestingly, we saw a relatively even distribution of up and down-regulated genes in bptf(RNAi) stem cells (Figs. 4A; 1603 up-regulated and 1886 down-regulated genes at day 9, 1288 up-regulated and 1332 down-regulated genes at day 12), which was somewhat unexpected given the overall loss of chromatin accessibility at gene promoters in these cells (Fig. 3E-F). However, it is possible that many of these DEGs are indirect effects of bptf(RNAi).

Fig. 4

Loss of BPTF leads to significant changes in gene expression that indicate the dysregulation of transcription and chromatin regulation. A) Volcano plots showing transcript changes in bptf(RNAi) stem cells versus control(RNAi) as detected by RNA-seq. Left plot shows data from stem cells isolated 9 days post-RNAi; right plot shows data from 12 days post-RNAi. Blue dots = significantly down-regulated (pAdj < 0.05; 1886 at day 9, 1332 at day 12), red dots = significantly up-regulated (pAdj < 0.05; 1603 at day 9, 1288 at day 12). Vertical dashed lines = log2FC +/-0.75. Horizontal dashed line = linear pAdj 0.05. B) Clustered heat map of day 9 and day 12 RNA-seq data in B; only genes with significant differential expression (pAdj < 0.05) at both time points are included. Changes in expression between time points correlate strongly (Supplemental Figure S3A). C) Venn diagrams comparing genes identified as differentially accessible (Fig. 3; csaw) and differentially expressed (A; DESeq2). Genes included have pAdj < 0.05. Overlap p-values were calculated using the Hypergeometric test. D) Plots comparing changes in chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) for those genes identified as significantly differential in both in C. Left plot shows correlation with RNA-seq data from day 9 post-bptf(RNAi) and right plot shows data from day 12. Correlation coefficients (R) and p-values were calculated using Spearman Rank correlation. The grey shaded areas around the regression lines represent 95% confidence intervals

We then asked if there was a correlation between BPTF-dependent changes in chromatin accessibility and gene expression. We first generated Venn diagrams to determine how many gene loci had both DA (Differentially Accessible) ATAC-seq peaks at their promoters and differential expression by RNA-seq (Fig. 4C). The overlap between DA promoters and differentially expressed genes (DEGs) was significantly greater than what is predicted by chance, as calculated with hypergeometric tests using the total number of annotated gene models as the universe. However, these comparisons do not assess the relationship between changes in chromatin accessibility and gene expression at individual loci. To address this question, we compared the log₂ Fold Change (log2FC) in chromatin accessibility (csaw) at gene promoters with their log2FC in gene expression (RNA-seq; Supplemental Figure S3C). When we included all genes with a calculated value in both datasets (i.e., any pAdj), we found no correlation between changes in chromatin accessibility and gene expression (Supplemental Figure S3C). When we restricted our analysis to genes with significant (pAdj < 0.05) changes in both datasets (Fig. 4D), we still found no significant correlation at day 9 post-RNAi and a marginally significant positive correlation at day 12. We also did not find a significant correlation when comparing differential expression with differential accessibility using peaks in other mapping categories (e.g., distal intergenic, 5’ UTR, etc.; data not shown). However, we note that those categories contain relatively fewer genes in each and it is more challenging to assign those peaks to the gene(s) they may be regulating.

Importantly, when we binned genes into four broad expression categories based on their average TPM (transcripts per million) in control(RNAi) cells and then compared the ATAC-seq signal at their promoters (+/- 1kB), we did observe a significant correlation between average gene expression and promoter accessibility (Supplemental Figure S3B). Moreover, RNAi of bptf had a significant impact on accessibility at genes of all expression levels except the highest bin (> 100 TPM). These data suggest that BPTF-dependent loss of chromatin accessibility can impact gene expression, but not at all gene loci and potentially not at those genes expressed at very high levels.

BPTF-mediated accessibility is enhanced at promoters with H3K4me3

Because BPTF is known to bind H3K4me3 and this modification is highly enriched at gene promoters [29, 88, 89], we wondered if its presence might distinguish those loci at which BPTF affects gene expression from those at which it does not. To address this question, we first examined H3K4me3 patterns in both control(RNAi) and bptf(RNAi) stem cells. We again isolated X1 stem cells from both control(RNAi) and bptf(RNAi) animals, isolated their chromatin, and performed ChIP-seq with an antibody specific to H3K4me3. As shown previously in S.mediterranea [29, 31] and many other organisms [88,89,90], nearly all MACS2-called H3K4me3 peaks were detected at gene promoters (Fig. 5A, B). Importantly, we did not observe a detectable change in H3K4me3 after knockdown of bptf (Fig. 5B), indicating that the chromatin accessibility changes we detected in bptf(RNAi) stem cells were not due to loss of its binding site. We then asked how many gene promoters have both a MACS2-called H3K4me3 peak and a MACS2-called ATAC-seq peak in control(RNAi) stem cells. We found that 28% of all promoters with a MACS2-called ATAC-seq peak +/- 1 kb also had an H3K4me3 peak within 1 kb (Fig. 5C). Conversely, 25% of all promoters with an H3K4me3 peak also had an ATAC-seq peak (Fig. 5C).

Fig. 5

H3K4me3 concentrates chromatin accessibility at promoters. A) Pie chart showing the distribution of mapped H3K4me3 peaks in control(RNAi) stem cells. B) Profile plot comparing H3K4me3 signal in control(RNAi) and bptf(RNAi) stem cells averaged across all genes. The shaded area represents the standard error. C) Venn diagram showing the number of promoters in planarian stem cells with an H3K4me3 peak, an ATAC-seq peak, or both. Overlap p-value was calculated using the hypergeometric test. D) Profile plot comparing H3K4me3 signal in control(RNAi) and bptf(RNAi) stem cells at the gene promoters in each group of the Venn in C (H3K4me3 only, H3K4me3 + ATAC peak, and ATAC peak only). The shaded area represents the standard error. E) Profile plot comparing ATAC-seq signal in control(RNAi) and bptf(RNAi) stem cells for gene promoters in each group of the Venn shown in C. The shaded area represents the standard error. F) Box plots measuring changes in ATAC-seq signal at the promoters of genes in each group in E. Statistical significance was determined using a Wilcoxon test (p-value with Bonferroni correction, **** = p ≤ 0.0001). G) Scatter plots comparing changes in chromatin accessibility (ATAC-seq; csaw) with changes in gene expression (RNA-seq; DESeq2) in bptf(RNAi) X1 stem cells at loci with an H3K4me3 peak at their promoters (left plot) versus those without an H3K4me3 peak (right plot). Data points shown have pAdj < 0.05 in both datasets. Correlation coefficients (R) and p-values were calculated using Spearman Rank correlation. The grey shaded areas around the regression lines represent 95% confidence intervals

The number of genes with both ATAC-seq and H3K4me3 peaks at their promoters may be fewer than expected based on other studies. However, we note the following: first, when we examined recently published H3K4me3 ChIP-seq and ATAC-seq data generated from whole planarian worms [31], it also showed a large number of ATAC-seq peaks at promoters that did not have H3K4me3 (Supplemental Figure S4C). Second, when we individually inspected several "ATAC-seq only" peaks using IGV, we found they had robust and specific signal near promoters (Supplemental Figure S5) and looked like true ATAC-seq peaks. Lastly, most comparable genomic datasets were generated from mammalian cell lines grown in culture conditions that were optimized for homogeneity [91, 92]. In contrast, our data was generated from isolated planarian stem cells that are known to be transcriptionally heterogenous and include a mix of pluri- and multipotent stem cells.

To determine whether the presence of an H3K4me3 peak impacts BPTF-mediated chromatin accessibility, we then examined the average H3K4me3 (Fig. 5D) and ATAC-seq signal (Fig. 5E) at genes annotated with: (1) an H3K4me3 peak only; (2) both an H3K4me3 peak and an ATAC-seq peak; and (3) an ATAC-seq peak only (groups based on Fig. 5C). This analysis revealed that bptf(RNAi) induces the clearest loss of accessibility at genes with both H3K4me3 and ATAC-seq peaks (Fig. 5E, middle), although it has a statistically significant impact on chromatin accessibility at loci in all three categories (Fig. 5F). Noticeably, genes with “ATAC-seq only” promoters had relatively higher levels of ATAC-seq signal overall (Fig. 5E, right) and the signal was highly variable between genes in this group (as indicated by the shaded standard error). We noted that this third group included significantly more “h1SMnG” gene models (i.e., putative non-coding) and wondered whether the chromatin state at these loci was different from that of “high confidence” gene models (defined in [31]). However, when we split the “ATAC-seq only” group into “high confidence” models versus “putative non-coding”, we did not see a major difference in the profile plot (Supplemental Figure S4D). We further noted that the average ATAC-seq signal at the “ATAC peaks only” loci was less concentrated at the TSS than that of gene loci with both ATAC and H3K4me3 peaks (Fig. 5E). Yet, as discussed above, the ATAC-seq only peaks did show robust signal to background (i.e., distinct peaks) when inspected at individual gene loci (Supplemental Figure S5). These data suggest that H3K4me3 may focus NuRF activity at more precise nucleosome positions relative to the TSS.

We then asked if genes with BPTF-dependent ATAC-seq peaks are more likely to exhibit changes in gene expression if their promoters are marked by H3K4me3 peaks. First, we used Venn diagrams to compare any significant change in accessibility (increased or decreased) with any significant change in gene expression (up or down) after bptf(RNAi) in X1 stem cells (Supplemental Figure S4F, G). Based on hypergeometric tests, changes in chromatin accessibility at any promoter (i.e., those with or without H3K4me3) have significant overlap with changes in expression. To determine if these correlations were also directional and at individual genes, we then performed correlation analyses. When we included all genes with a calculated log2FC (any pAdj) for both ATAC-seq and RNA-seq, neither genes with H3K4me3 at their promoters nor those without showed a statistically significant correlation (Supplemental Figure S4H). However, when we restricted our analysis to genes with significant changes (pAdj < 0.05) in both ATAC-seq and RNA-seq after bptf(RNAi), there was only a significant correlation for those genes with H3K4me3 peaks at their promoters but not those without (Fig. 5G).

Set1 targets are significantly regulated by BPTF/NuRF

We then wondered if the type of H3K4me3 peak at a gene promoter impacted the role of BPTF. To address this question, we first realigned our previously published H3K4me3 ChIP-seq data generated from set1(RNAi) and mll1/2(RNAi) X1 stem cells on the new S.mediterranea genome assembly (Supplemental Figure S6 and Supplemental Table 3). After identifying the new coordinates of Set1 and MLL1/2-dependent target loci, we then compared the ATAC-seq signal at these distinct sets of gene loci (Fig. 6A). Unsurprisingly given our previous findings [29], Set1 target genes had considerably more ATAC-seq signal than MLL1/2 target genes (Fig. 6A, B). At the same time, both Set1 and MLL1/2 target genes lost a statistically significant amount of ATAC-seq signal in bptf(RNAi) stem cells (Fig. 6B), even though the chromatin at MLL1/2 target genes was overall less accessible.

Fig. 6

BPTF shares functional genomic targets with Set1. A) Profile plot comparing ATAC-seq data from control(RNAi) and bptf(RNAi) stem cells at Set1 (left) and MLL1/2 (right) gene targets. The shaded area represents the standard error. B) Quantitation and statistical analysis of ATAC-seq data in A; statistical significance was determined using a Wilcoxon test (p-value with Bonferroni correction). * = p ≤ 0.05, ** = p ≤ 0.01, *** = p ≤ 0.001, **** = p ≤ 0.0001). C) Venn diagram comparing genes identified as likely targets of Set1, MLL1/2, and BPTF in X1 stem cells. Targets = genes with loss of H3K4me3 or ATAC-seq at promoters (+/- 1 kb of TSS) in their respective RNAi conditions. p-val = 8.65e-74 for BPTF targets overlap with Set1 targets; p-val = 1.04e-75 for BPTF targets overlap with MLL1/2. D) Venn diagrams comparing common Set1 and BPTF targets (217 genes) with DEGs at day 12 (left) as well as common MLL1/2 and BPTF targets with DEGs at day 12 (right). P-values were calculated using a hypergeometric test. E) Plots comparing differential changes in chromatin accessibility (ATAC-seq) with differential changes in gene expression (RNA-seq) in bptf(RNAi) stem cells at Set1 versus MLL1/2 target genes (from C, D). Data points shown have pAdj < 0.05 in both datasets. Correlation coefficients (R) and p-values were calculated using Spearman Rank correlation. The grey shaded areas around the regression lines represent 95% confidence intervals. F) GO Term enrichment analysis for genes with reduced chromatin accessibility in bptf(RNAi) X1 stem cells (“BPTF-dependent loci”), reduced H3K4me3 in set1(RNAi) X1 stem cells (“Set1-dependent loci”), and reduced H3K4me3 in mll1/2(RNAi) X1 stem cells (“MLL1/2-dependent loci”)

We then asked if bptf(RNAi) had different impacts on gene expression at Set1 versus MLL1/2 target genes. First, we compared promoters at which the chromatin state was affected by depletion of either Set1, MLL1/2, or BPTF (Fig. 6C). As expected, Set1 and MLL1/2 target distinct gene loci and a subset of each are affected by loss of BPTF. We then compared all common Set1 + BPTF target genes (217 genes) with genes differentially expressed in bptf(RNAi) stem cells (Fig. 6D, left). The overlap of these gene sets was significantly greater than what is predicted by chance (hypergeometric test). This was also true for the parallel comparison of genes targeted by both MLL1/2 and BPTF (57 genes; Fig. 6D, right). However, these analyses did not consider the degree or direction of change (i.e., increased or decreased accessibility or expression). When we compared significant ATAC-seq changes with significant RNA-seq changes at individual genes, only genes targeted by both Set1 and BPTF (but not MLL1/2 + BPTF) showed a significant, positive correlation (Fig. 6E). Notably, this was only true for RNA-seq changes at day 12 post RNAi, suggesting the effects of changing chromatin accessibility do not immediately impact stable transcript levels. It was also only true if we restricted the analysis to significant values (i.e., pAdj < 0.05 in both RNA-seq and ATAC-seq datasets; Fig. 6E) but not if we included log2FC values with any pAdj (Supplemental Figure S6C, D). Separately, when we performed Gene Ontology (GO) term analysis for genes targeted by BPTF, Set1 and MLL1/2 (Fig. 6F), we found that several molecular function terms were enriched in both BPTF and Set1 lists (e.g., DNA binding, structural constituent of chromatin). In contrast, terms enriched in the descriptions of MLL1/2 target genes were largely related to cilia.

We previously showed that Set1 and MLL1/2-dependent loci are marked by measurably different H3K4me3 patterns; Set1 catalyzes peaks with significantly wider widths than those added by MLL1/2 [29, 86] (Supplemental Figure S6A, B). To test if BPTF/NuRF activity is affected by different H3K4me3 peak widths, irrespective of the specific KMTase creating them, we divided all H3K4me3 peaks from control(RNAi) stem cells into three categories: “narrow”, “medium”, and “broad” (Supplemental Figure S7A-C). After plotting the average ATAC-seq signal in each group (Supplemental Figure S7D), we found that the overall ATAC-seq signal increased with H3K4me3 peak width (Supplemental Figure S7E, F). Unexpectedly, when we analyzed the reduction of BPTF-dependent chromatin accessibility across these three groups, the loss of ATAC-seq signal at genes with broad H3K4me3 peaks was not statistically significant (Supplemental Figure S7E, F). This was unexpected, as the average width of these H3K4me3 peaks was similar to those marked by Set1 (Supplemental Figure S6B, S7B). These findings suggest that Set1 and MLL1/2 target loci have additional features, beyond the H3K4me3 peaks they catalyze, that recruit, enhance and/or restrict BPTF/NuRF function (Fig. 7).

Fig. 7

Model of BPTF/NuRF function in planarian stem cells. A) SMED-BPTF is predicted to bind histone H3 lysine 4 trimethylation (H3K4me3; orange circles), recruiting the rest of the NuRF complex and its ATP-dependent chromatin remodeling activity to transcription start sites. Many other genes in the planarian genome are not marked by H3K4me3, but do have BPTF-dependent regions of chromatin accessibility. It is possible that other chromatin modifications, such as acetylation (magenta triangles), are responsible for recruiting and/or stabilizing BPTF/NuRF. B) Genes marked by Set1-dependent H3K4me3 have more open chromatin and higher expression than those marked by MLL1/2-dependent H3K4me3

Discussion

Planarian stem cells have many fascinating cellular characteristics and are required for two extraordinary features of their organism: (1) planarians can repeatedly and reproducibly regenerate complex tissues and organs, and (2) they are effectively immortal. Although many genes have been identified as essential regulators of planarian stem cell function, including genes that encode chromatin proteins [21,22,23,24, 26, 28, 29], the molecular mechanisms and genetic targets through which they accomplish this regulation are not well understood. Here, we examined the role of the chromatin remodeling complex NuRF in planarian stem cell biology. By knocking down its largest and specific subunit, BPTF, we both targeted the NuRF component that provides its chromatin-state specificity [37, 93] and tested its effects on transcription in vivo. Our findings show that BPTF, likely through the remodeling activity it recruits to specific gene loci, is critical for planarian stem cell function. We are currently working on analyses to examine the homology and function of these BPTF-sensitive genes, many of which encode proteins with identifiable domains (e.g., zinc fingers) but unclear homology.

Limitations of the study

When interpreting our data, several points must be considered. First, we acknowledge that the method for isolating stem cells from planarians is based on DNA content, i.e., cells with > 2n DNA [59]. This means that most cells in our study, and in most molecular/genomic studies of planarian stem cells, are in the late S/G2/M phases. It is possible that this characteristic may skew our chromatin accessibility data (ATAC-seq) toward a less accessible state, as chromosomes become highly compacted in M phase and studies using synchronized cells have reported that ATAC-seq signal decreases as cells progress through these phases of cell cycle [94]. In this study, we are comparing cells in which we have depleted an integral subunit of a chromatin remodeling complex, BPTF, to control cells. Despite the possible dampening of ATAC-seq signal due to cell cycle-induced compaction, we still detect many significant changes in chromatin accessibility after its depletion. However, it is possible that we would detect many more if we could isolate all stem cells, including those with 2n DNA content (in G1 or G0).

Second, knockdown of bptf may cause a small, but consistent, loss of accessibility at many gene promoters, but these changes may not reach the threshold required by csaw to identify them as differentially accessible (DA) peaks. Related, some gene loci may have highly dynamic chromatin states, making it difficult to assay their accessibility reproducibly across replicates. Both scenarios could depress our ability to detect significant changes in chromatin accessibility in bptf(RNAi) stem cells, which then also affects our correlation analyses with gene expression. Third, we are using standard RNA-sequencing measurements, which measure all mRNA transcripts in our cells. As these methods largely detect steady-state cytoplasmic mRNA transcripts, they are not an ideal readout of transcriptional activity. Assays such as GRO-seq or NAC-seq would be better measurements of dynamic transcriptional changes, but they are not currently tractable in the planarian model. Nonetheless, these issues are not unique to planarian studies; similarly weak correlations between changes in ATAC-seq and RNA-seq have been reported in other species and contexts, suggesting they may suffer from similar limitations [95].

Finally, the planarian field is still in the early stages of understanding the genomes of these animals and how they are regulated. Chromosome-level genome assemblies for S.mediterranea were only recently released [31, 96], and work remains to be done to understand its regulatory networks. For example, although the S.mediterranea genome is relatively large (~ 800 MB) with sizeable regions of repetitive sequence, the genic regions are often very dense. These details are relevant to any genomic assay, as they impact steps like mapping ATAC-seq peaks to their nearest gene model. Nevertheless, we anticipate that our data will provide a meaningful contribution to the characterization of gene loci that regulate planarian stem cell function.

Conclusions

Using ATAC-seq, we identified many specific regions of the planarian genome at which chromatin accessibility changed significantly after bptf knockdown. We confirmed that, as expected, the chromatin at Set1 target genes is significantly more accessible than that at MLL1/2 targets. Moreover, genes with both a BPTF-dependent loss of accessibility and a Set1-dependent H3K4me3 peak showed the strongest overall correlation with transcription. Finally, we observed that knockdown of bptf replicates many of the physiological, cellular, and molecular changes seen in set1(RNAi) planarians, including the inability of their stem cells to repopulate after treatment with ionizing radiation.

Our study is also important because it likely informs how BPTF operates in vivo. In humans, BPTF appears to have dose-dependent functions: haploinsufficiency of BPTF causes abnormal craniofacial and brain developmental [97, 98] whereas over-expression of BPTF is linked to oncogenic behavior in various cancers [99,100,101]. These data suggest that some cell types and/or genomic loci are more sensitive to the amount of BPTF-dependent NuRF activity than others. However, the exact mechanisms driving these functional outcomes in different in vivo settings remain unclear. Because we robustly reduced, but did not eliminate, bptf expression, we identified the loci and genes that are most sensitive to precise levels of NuRF activity in planarian stem cells. Despite this incomplete loss, bptf(RNAi) animals have severe and penetrant phenotypes, showing that normal NuRF activity levels are critical for in vivo cellular function in planarians as in humans. Future studies will be aimed at understanding the functions and gene regulatory networks of those most sensitive target genes.