Drosophila Brahma complex remodels nucleosome organizations in multiple aspects

ATP-dependent chromatin remodeling complexes regulate nucleosome organizations. In Drosophila, gene Brm encodes the core Brahma complex, the ATPase subunit of SWI/SNF class of chromatin remodelers. Its role in modulating the nucleosome landscape in vivo is unclear. In this study, we knocked down Brm in Drosophila third instar larvae to explore the changes in nucleosome profiles and global gene transcription. The results show that Brm knockdown leads to nucleosome occupancy changes throughout the entire genome with a bias in occupancy decrease. In contrast, the knockdown has limited impacts on nucleosome position shift. The knockdown also alters another important physical property of nucleosome positioning, fuzziness. Nucleosome position shift, gain or loss and fuzziness changes are all enriched in promoter regions. Nucleosome arrays around the 5′ ends of genes are reorganized in five patterns as a result of Brm knockdown. Intriguingly, the concomitant changes in the genes adjacent to the Brahma-dependent remodeling regions have important roles in development and morphogenesis. Further analyses reveal abundance of AT-rich motifs for transcription factors in the remodeling regions.


INTRODUCTION
Eukaryotic genomic DNA is densely packaged in nucleus through high-order chromatin architecture. Nucleosomes are the fundamental structural unit of chromatin, composed of a histone octamer wrapped by ∼147 bp DNA (1). Nucleosome positioning plays a critical role in numerous biological processes mainly through regulating the accessibility of DNA. There are various factors affecting nucleosome positioning such as DNA sequence, histone variants, histone tail modifications, ATP-dependent chromatin remodeling enzymes and histone chaperones. It has been reported that the AA/TT dinucleotides occurred in a biased and/or periodic arrangement across nucleosomal DNA around the 5 end of genes in yeast (2). In contrast, the nucleosome positioning sequence pattern in Drosophila is CC/GG (3). Both histone variants and histone modifications often exist in the nucleosomes around the 5 end of genes that could facilitate nucleosome eviction (4).
The SWItch/Sucrose NonFermentable (SWI/SNF) family, one of the ATP-dependent chromatin remodeling complex families, has diverse roles. Conditional knockout study in mice showed that Brg1, a member of the SWI/SNF family, played an important role in neural development through regulating sonic hedgehog signaling (5). The ablation of Brg1 in neural stem cells in mice failed to maintain neural stem cells in a gliogenic state and switched from neurogenesis to gliogenesis (6). Other studies implicated that SWI/SNF complexes also functioned in self-renewal and pluripotency of mouse embryonic stem cells, heart and thymocyte development (7)(8)(9). Inactivation of the key subunits of the mouse SWI/SNF complex (Snf5 or Brg1) resulted in a loss of nucleosome occupancy at target promoters (10). In Drosophila, there are two distinct SWI/SNF complexes, Brahma-associated proteins (BAP) and polybromocontaining BAP (PBAP), that are characterized by specific accessory subunits. Studies have found that the Drosophila Brahma complex played a role in wing development (11), cell cycle control (12), self-renewal of neural stem cells (13) and intestinal stem cell proliferation (14). The signature subunit Osa of BAP complex was required to repress the expression of Wingless target genes (15). The subunit Osa was also needed for the Drosophila wing development by regulating the epidermal growth factor receptor (EGFR) pathway (16) and mediating the expression of Apterous target genes (17). The signature subunit Bap170 of PBAP complex was required in normal eggshell development (18). Snr1, a highly conserved subunit of Brahma complex also functioned in the development of wing vein (19).
In spite of recent substantial effort, the role of Brahma complex, the common ATP-dependent catalytic subunit of Drosophila SWI/SNF class of chromatin remodelers, in affecting nucleosome landscape in vivo is unclear. To gain insights into this question, we employed the UAS-Gal4 system driven by the promoter of heatshock protein in Drosophila. We optimized a heating program to precisely control the timing of Brm knockdown in Drosophila third instar larvae. Then we generated and compared the highresolution maps of nucleosome positions in the genome before and after Brm knockdown. Our results showed that knockdown of Brm led to extensive changes in nucleosome landscape. The resultant remodeled regions were enriched in the promoter regions. The differentially expressed genes adjacent to the remodeled regions had functions for disc or tissue morphogenesis.

Fly stocks and crossing
Stocks used in this paper are UAS-mCD8GFP, UAS-Brm-IR (#31712) and Hs-Gal4 from Bloomington Stock Center (http://flystocks.bio.indiana.edu). All transgenes are on the third chromosome. We crossed UAS-Brm-IR males with Hs-Gal4 females to generate the transgenic strain containing UAS-Brm-IR and Hs-Gal4 for Brm knockdown (experimental group, short for BrmIR). Similarly, the transgenic strain containing UAS-mCD8GFP and Hs-Gal4 (control group, short for GFP) was generated by crossing UAS-mCD8GFP males with Hs-Gal4 females and expected to remove the influence of UAS-Gal4 system.

RNAi by heatshock at the third larvae stage of D. melanogaster
We collected the fly fertilized eggs of the control and experimental group, respectively, and cultured them at normal temperature (25 • C) for 83.5 h (i.e. the third larvae stage). Next, we cultured the larvae at 37 • C for 0.5 h. During the heatshock, Brm was knocked down by the expression of Brm inverted repeat (IR) transgene in the experimental group (denoted as BrmIR37). In contrast, Brm was not knocked down in the control group (denoted as GFP37) because of no expression of Brm IR transgene. Then we recovered the larvae at 25 • C for 4.5 h. Repeat the heatshockand-recovery treatment twice, but with recovery time of 3.5 and 2 h, respectively. After the heatshock treatment, the larvae were collected and used in this study. It is noteworthy that we tried other heatshock programs, e.g. heatshock in an earlier stage or longer heatshock time, leading to embryonic development arrest at early stages. The current heatshock program significantly reduced the Brm transcription and the embryos could develop until the late pupae. Therefore, this heatshock program was used in our study.

Total RNA extraction and qRT-PCR
Total RNA was extracted from the larvae using TRIzol (Invitrogen) that underwent the above heatshock treatment. Reverse transcription was performed on 500 ng of total RNA in 10 l reaction system for each sample using PrimeScript RT Master Mix Perfect Real Time Kit (TAKARA, Code: DRR036A). The realtime polymerase chain reaction (PCR) was carried out in triplicate using F-416L DyNAmoTM ColorFlash SYBR R Green qPCR Kit (Thermo). Relative quantitative analysis of Brm expression levels was calculated using on 2 − Ct formula described in a previous study (20) with gene tub as the reference gene. Primers for Brm are sense 5 -CCAATGCCGAGCGTGAAC-3 and antisense 5 -ACTCATCTGTCTGCGACAGTAGGA-3 .

RNA-seq data analysis
The RNA sequencing libraries were constructed from the extracted RNA using standard Illumina libraries prep protocols. RNA-seq was performed on Illumina HiSeq2000 platform. Sequencing reads were aligned to annotated Drosophila transcripts (FlyBase r5.43) using TopHat with zero mismatch. The statistics of raw reads, mapped reads were summarized in Supplementary Table S1. The differentially expressed genes were identified by the tool Cuffdiff. For each transcript, reads per kilobase per million (RPKM) mapped reads were calculated to evaluate the expression level.
The pellet was washed with 500 l MNase digestion buffer (10 mM Tris-HCl (pH 7.5), 15 mM NaCl, 60 mM KCl, 1 mM CaCl 2 , 0.15 mM spermine, 0.5 mM spermidine, 1×EDTA-free protease inhibitor cocktail) and centrifuged at 16 000g for 10 min at 4 • C. The nuclei were resuspended in 500 l MNase digestion buffer plus 15 U MNase (Worthington, LS004797) at 37 • C for 30 min. The MNase digestion was terminated on ice by adding EDTA to a final concentration of 10 mM for 10 min. The mixture was centrifuged at 16 000g for 10 min at 4 • C. The supernatant was discarded. The pellet was resuspended in 500 l A3 buffer (140 mM NaCl, 15 mM HEPES (pH 7.6), 1 mM EDTA, 0.5 mM EGTA, 1% Triton X-100, 0.1% sodium deoxycholate, 0.1% SDS, 1×EDTA-free protease inhibitor cocktail). RNA was removed by digestion with 5 l 10 mg/ml RNaseA at 37 • C for 1 h in waterbath and with 3 l 20 mg/ml Proteinase K at 65 • C overnight in waterbath, respectively. Nucleosomal DNA was extracted by phenol-chloroform and dissolved in ddH 2 O. Mononucleosomal DNA fragments were extracted on a 2% agarose gel (Supplementary Figure S1) and purified by the MinElute Gel Extraction Kit (QIAGEN, 28604). The purified mononucleosomal DNA was subjected to massively parallel DNA sequencing on Illumina HiSeq2000 using single end protocol.

Nucleosome prediction and analysis of positioning dynamics
Sequencing reads were aligned to Drosophila melanogaster reference genome (dm3) using Bowtie with up to two mismatches. The statistics of raw reads, mapped reads were summarized in Supplementary Table S1. Only the uniquely mapped reads were retained in this study. The resulting sequence read distribution was used to identify nucleosomes using the peak-calling tool GeneTrack (22,23). Each nucleosome was assigned to either of promoter, genic or intergenic regions depending on in which region the midpoint of the nucleosome located. The nucleosome coordinates were given in Supplementary data file 1 (GFP37) and file 2 (Br-mIR37).
We located the closest nucleosome in GFP37 sample for each nucleosome in BrmIR37 sample. The two nucleosomes were defined the same nucleosome if the distance between their midpoint was <80 bp. Otherwise, nucleosome gain and loss occurred. For each pair of nucleosomes, position shift was calculated by subtracting the midpoint coordinate of the nucleosome in GFP37 from that in BrmIR37. Nucleosome fuzziness was calculated as the standard deviation of the midpoint of the set of reads defining the same nucleosome as described previously (3). It measures how spread out a nucleosome position is. Nucleosome fuzziness difference was calculated by subtracting the nucleosome fuzziness in GFP37 from that in BrmIR37.

Nucleosome occupancy change analysis
We scanned the Drosophila genome with a 200-bp window, and calculated nucleosomal sequencing read count (RPKM) in each window as its nucleosome occupancy. The ratio of nucleosome occupancy after to before Brm knockdown in each window represents the nucleosome occupancy change. We retained the windows with at least 2-fold nucleosome occupancy change, i.e. ratio ≥2 (increase) or ratio ≤0.5 (decrease), for the next analysis. Each window was assigned to either of promoter, genic or intergenic region depending on which genomic region overlaps with ≥50% of the window length. Then, the enrichment of nucleosome occupancy increase in each genomic region is the total number of nucleosome occupancy increase windows normalized by the total length of corresponding genomic regions. The enrichment of nucleosome occupancy decrease was analyzed in the same way.

Nucleosome organization change around TSS
The annotation of all Drosophila genomic features were downloaded from FlyBase release 5.43 (ftp://ftp.flybase.net/ genomes/dmel/dmel r5.43 FB2012 01/fasta/). The nucleosome array within ±1 kb of the Transcription Start Site (TSS) was plotted as the previous study (24). Briefly, nucleosomes within ±1 kb of the TSS formed the nucleosome array of a gene. The nucleosome position was defined by its dyad. The length was equal to the fuzziness value that measured the spread-out of a nucleosome. We compared the nucleosome array of each gene before and after Brm knockdown. As a result, each site had either of the two possible nucleosomal states: occupied by or depleted of a nucleosome. There was no difference in nucleosome organization at a site that had the same nucleosomal state before and after Brm knockdown. Otherwise, nucleosomal state at a site was changed after Brm knockdown. Then nucleosome organization difference was clustered by K-means (K = 5). Note that we tried several K values and K = 5 gave distinct patterns among clusters and least variation within each cluster. The clustering results were plotted in heatmap.

Composite distribution of nucleosome relative to TSS
The composite nucleosome occupancy landscapes were calculated by aggregating read count at each distance relative to the TSS as described previously (3). Briefly, given a subset of TSSs, GFP37 and BrmIR37 nucleosome reads within ±1 kb of each subset of TSSs were collected, respectively. Each read represented a nucleosome. Thus, we extended each read toward 3 end to a length of 147 bp. The midpoint of extended read defined the nucleosome position. Nucleosome distances to the TSS were binned in 10-bp intervals. Nucleosome occupancy was the total read count in the corresponding bin. Finally, nucleosome occupancy in each bin was further normalized as RPKM, and smoothed with five bins.

Brahma-dependent chromatin remodeling regions
Global Brahma-dependent chromatin remodeling regions were detected by the similar method used in the previous study (25). In brief, the normalized nucleosome signal of each genomic site, N i , is the ChIP tag count at coordinate i normalized by total tag count in the sample. Then the nucleosome signal change of site i, S i , was measured by the log 2 transformed ratio of N i (BrmIR37) to N i(GFP37) . The average length of linkers in Drosophila is 18 bp (3). Thus, the midpoint distance of two adjacent nucleosomes is 165 bp on average. To reduce the false positive rate, we required a chromatin remodeling region that span at least three nucleosomes, and defined remodeling score at coordinate i as Z i = S i (S i−165 + S i+165 ). The sites whose absolute value of the Z score is ≥1 (FDR = 0.05) were defined as remodeling blocks. Blocks <165 bp apart were merged to one remodeling region. This analysis identified 3758 distinct regions, typically of 500-700 bp in length.

Accession numbers
The RNA-seq and MNase-seq data sets have been deposited in ArrayExpress database under accession numbers E-MTAB-1966 and E-MTAB-1967, respectively.

Brm knockdown leads to perturbation on gene transcription with diverse functions
We generated the control fly strain (UAS-mCD8GFP/Hs-Gal4) and the Brm RNAi fly strain (UAS-Brm-IR/Hs-Gal4). The Brm was knocked down by the expression of Brm IR transgene through heatshock (see 'Materials and Methods' section for the details). The qRT-PCR results showed that the expression level of Brm was significantly decreased by knockdown ( Figure 1A). The western blotting result further confirmed the largely reduced level of Brahma complex in the knockdown sample ( Figure 1B). The normalized read count on Brm gene body confirmed the high expression level of Brm IR transgene and low expression level of Brm in the Brm RNAi fly strain when compared to the control fly strain ( Figure 1C). This suggested that the UAS-Gal4 system worked well to knockdown Brm in vivo in our study. To investigate the effect of Brm knockdown on global gene expression, we identified 872 significantly differentially expressed genes whose expression change was >2-fold after Brm knockdown ( Figure 1D). The differentially expressed genes included the hormone-responsive Ecdysone-induced genes (Eig) and homeotic genes that were reported to be the targets of Brahma in the previous studies (26)(27)(28). Gene Ontology (GO) analysis found that the differentially expressed genes had diverse functions including molting cycles, de novo protein folding, immune response, sensory perception and more ( Supplementary Figure S2). A previous study showed that the Brm mRNA level was significantly higher in 0-12 h Drosophila embryos than at later developmental stages (27). Consistent with it, knockdown of Brm in early embryos arrested the development of Drosophila embryos in our study. Moreover, Brm depletion from the zygote was lethal at late stages of embry-onic development (29). Therefore, Brm has vital functions for the development of Drosophila embryos.

Brahma alters global nucleosome organization and the physical properties associated with nucleosome positions
We aligned the nucleosomal reads from MNase-seq against Drosophila genome. The read count in a given region indicates its nucleosome occupancy. We compared the genomewide nucleosome occupancy before and after Brm knockdown. The results showed that nucleosome occupancy changes spread out the entire genome after Brm knockdown (Figure 2A and Supplementary Figure S3). To examine whether there is a bias in the distribution of nucleosome occupancy change in different genomic regions, we calculated the enrichment of nucleosome occupancy change in promoter, genic and intergenic regions. Our results showed a similar level of nucleosome occupancy increase in the all three genomic regions. In contrast, promoter and intergenic regions had higher enrichment of nucleosome occupancy decrease than genic regions. Moreover, nucleosome occupancy decrease was sufficiently higher than occupancy increase in either of the genomic regions ( Figure 2B). This was partially attributed to the fact that Brahma-containing remodeling complexes (P)BAP promoted formation of nucleosome at their target sequences (30).
To further examine the impact of Brm knockdown on the physical properties associated with nucleosome positions, we employed GeneTrack (23) to generate genomewide maps of mononucleosome positions before and after Brm knockdown for comparison, respectively. We first investigated how Brm knockdown affected nucleosome position shift. 89.6% of total nucleosomes (406 286) shifted after Brm knockdown. More specifically, 66.4% of nucleosomes shifted <10 bp. There was no shift in only 8.5% of nucleosomes ( Figure 2C). Thus, Brm knockdown led to a limited change in nucleosome positions. Unlike occupancy change, these shifted nucleosomes were enriched in promoter regions and gene body ( Figure 2D). In addition to nucleosome shift, Brm knockdown also caused nucleosome disassembly (loss) and assembly (gain). Intriguingly, there was also an enrichment of nucleosome gain or loss in promoter regions and gene body ( Figure 2E). There was similar number of nucleosome gain and loss within each genomic regions. Although nucleosome gain or loss occurred in 7381 genes, only 379 (5%) genes were significantly differentially transcribed. These findings suggested that nucleosome organization changes was not tightly correlated with gene transcription changes. The similar findings were also reported in yeast (31). GO term analysis found that these genes had a key role in positive regulation of kinase activity (Supplementary Figure S4).
Fuzziness is another important physical property of nucleosome positions. It measures the delocalization of a nucleosome position. The larger the fuzziness, the more delocalized a nucleosome positioning is. Conversely, a highly phased nucleosome has a very small fuzziness value and its position is fixed irrespective of cell state or developmental stage. Interestingly, Brm knockdown significantly changed the fuzziness of 19 081 nucleosomes (Student's t-test, P-value < 0.01). Consistent with nucleosome shift and gain/loss, nucleosomes with fuzziness change were enriched in promoter and gene body ( Figure 2F). The fuzziness change was bidirectional. Some nucleosome positioning became delocalized after Brm knockdown and some was opposite. Unlike nucleosome occupancy change, there was no bias in fuzziness increase or decrease.

Nucleosome landscape changes around TSS regions arise from knockdown of Brm
The −1, NFR (nucleosome free region), +1, +2, +3, etc. canonical nucleosome arrangement around TSS regions plays an important role in gene transcription regulation and is linked to many biological processes (3,4,22,32,33). To gain insights on how Brm knockdown alters nucleosome profiles around TSS regions, we explored nucleosome organization difference in the regions surrounding TSS. The unsupervised clustering of nucleosome organization changes around TSS regions revealed five distinct patterns (Figure 3). The original composite distribution of nucleosomes around TSS of approximately two-thirds of genes (Cluster I) showed a very similar pattern before and after Brm knockdown. There was only slight shift in +3, +4, +5, etc. downstream nucleosomes and the nucleosomes >500 bp upstream of TSS (Cluster I). In contrast, Brm knockdown resulted in pronounced shift to 3 end in +3, +4, +5, etc. downstream nucleosomes. The positioning phase of highly phased −1, +1 and +2 nucleosomes remained un-changed (Cluster III). Another pattern of nucleosome profiles showed position shift of nucleosomes nearby 1-kb upstream of TSS (Cluster VI). There was also marked shift of the entire canonical nucleosomes around TSS (−1, +1, +2, +3, +4, +5, etc. nucleosomes). The shift direction was towards 5 end in a group of genes (Cluster II) and toward 3 end in the other group of genes (Cluster V). The functional analysis showed that the different classes of genes had the enrichment of distinct GO terms ( Supplementary Figure S5). For example, respiratory system development, cell motion, cuticle development, etc. were enriched in Cluster II genes. Oxidation reduction, mesoderm development, gastrulation, etc. were enriched in Cluster V genes.

AT-rich motifs for transcription factors are abundant in Brahma-dependent chromatin remodeling regions
To systematically identify Brahma-dependent chromatin remodeling regions, we developed a comparative approach to detect difference in nucleosome occupancy before and after Brm knockdown. The analysis identified 3758 remodeling regions that were enriched in the promoter regions. Particularly, there was abundance of the remodeling regions adjacently upstream of TSS ( Figure 4A and B). We further collected 1301 genes whose TSS located within 1 kb of these remodeling regions. Only 14.3% of these genes were differentially expressed. Seventy genes were up-regulated and 116 genes were down-regulated after Brm knockdown ( 4C). Intriguingly, changes of these gene expression affected embryonic development and morphogenesis ( Figure 4D).
The key function of nucleosome positioning is to occlude many regulatory DNA elements. Therefore, we scanned the remodeling regions to identify motifs as the potential binding sites for transcription factors (TFs) using Cistrome (34). The motif discovery results found a bunch of motifs for TFs (Supplementary Table S2). This suggested that Brahma complex regulate the accessibility of these motifs through chromatin remodeling. Interestingly, sequence composition analysis revealed that the sequences of these motifs contained an extremely high level of A and T nucleotides (Figure 4E). The previous work have reported that the rigid poly (dA:dT) tracts disfavored nucleosome formation (35,36). Coincidently, the target sites of Drosophila SWI/SNF com-plexes disfavored nucleosome formation (30). This partially explained Brm knockdown resulted in more nucleosome occupancy decrease than nucleosome occupancy increase (Figure 2A and B). Notably, presence of TF motifs in the remodeling regions provided a mechanism by which TFs helped recruit chromatin remodelers to their target sites. For example, yeast ISW2 chromatin remodeling complex repressed early meiotic genes by establishing an inaccessible chromatin structure upon recruitment by Ume6p (37).

DISCUSSION
Nucleosomes are the primary repeating units of eukaryotic chromatin structure, controlling access of DNA. Thereby, nucleosome organization plays important roles in many biological processes. ATP-dependent chromatin remodeling Chromatin remodelers counteract DNA sequence-driven nucleosome positioning through binding to their target sites. Interestingly, ChIP-chip analysis found the binding overlapped between chromatin remodeling complexes (P)BAP, NURD and INO80, but not ISWI (30). We also compared our Brahma-dependent chromatin remodeling regions with the reported binding sites of (P)BAP, NURD, INO80 and ISWI (30). Intriguingly, the remodeling regions were enriched on the binding sites of (P)BAP, NURD and INO80, but depleted on the binding sites of ISWI (Supplementary Figure S6). This results could be partially attributed to the fact that the target sits of (P)BAP, NURD and INO80 disfavored nucleosome position whereas the remodeling sites of ISWI promoted nucleosome placement (30). The enrichment of Brahma-dependent chromatin remodeling regions on the binding sites of (P)BAP, NURD and INO80 may also implicate their potential involvement in Brahma-dependent chromatin remodeling. However, more validation assays are required to resolve whether there is an interplay between the different families of chromatin remodelers during Brahma-dependent chromatin remodeling. For example, the identification of the binding sites of Brahma by ChIP-seq in larvae and overlapping analysis of their binding sites.
Chromatin remodeling changes nucleosome positioning and consequently gene expression. However, the extent of impact on gene expression of nucleosome positioning de- pends on the genomic site where nucleosome remodeling occurs. It was reported that the canonical arrangement of −1, NFR, +1, +2, +3, etc. nucleosomes around TSS played important roles in gene expression (3,4). The active genes possessed this canonical nucleosome organization around their TSSs whereas the TSSs of silent genes were occupied by a nucleosome blocking their access (33,38). In our study, 14.3% of the genes whose TSS located within 1 kb of the remodeling regions expressed differentially after brm knockdown whereas only 5% genes with nucleosome gain or loss on the gene body expressed differentially. There was significantly over proportion of the genes nearby the remodeling regions that expressed differentially after brm knockdown (χ 2 test, P-value = 2.45E−29). This further con-firmed that nucleosome occupancy change in TSS had a stronger correlation with gene expression. Although nucleosome occupancy usually occludes access to DNA sequences, regulatory DNA elements on nucleosome surface could also be bound by TFs through nucleosomal rotational setting (22). Thus, the occlusion effect of nucleosome occupancy on DNA access was weakened. Of note, in addition to nucleosome occupancy, other factors including histone modifications and DNA methylation also regulate gene expression. As a matter of fact, these factors cooperated to regulate gene expression through complex cross talks (39). Therefore, there was not always a strong association between nucleosome remodeling alone and gene expression change, and vice versa. For example, inactivation of Isw1 and Chd1 in yeast obtained the similar findings that the resulting perturbation in nucleosome organizations was not correlated with changes in gene transcription (31). In summary, Brahma-dependent chromatin remodeling alters nucleosome landscape and regulates gene transcription in an intricate way.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.