PARP1 is involved in mRNA splicing
We previously showed that PARP1 KD in Drosophila S2 cells results in changes in alternative splicing of several genes [5]. Our goal in this study is to understand mechanistically how PARP1 modulates chromatin structure to regulate splicing decisions. We chose to analyze this mechanism in depth at two genes—AKAP200 (hence called AKAP) and CAPER—because we found that (1) PARP1 binds to nucleosomes within their gene bodies and (2) PARP1 depletion correlated with changes in the splicing decisions observed through RNA-seq genome-wide analyses [5] (Fig. 1a). We used RNA interference to deplete PARP1 in S2 Drosophila cell lines (Additional file 1: Fig. S1), to assess whether PARP1 binding on AKAP and CAPER exon–introns reflects a direct role for PARP1 in alternative splicing decisions at these genes. We then performed PCR with exon junction spanning primers to validate the splicing changes (Fig. 1b). These results validated the genome-wide studies [5] and showed that in the absence of PARP1, differentially spliced transcripts for AKAP and CAPER were produced. Furthermore, a second siRNA (siRNA2) targeting a different region of PARP1 (Additional file 1: Fig. S1A and B) confirmed these results (Additional file 1: Fig. S2). To test whether this effect is due to PARP1 directly or its enzymatic activity, we inhibited its PARylation using PJ34 (Additional file 1: Fig. S1C) and showed that PARylation inhibition effected no changes in splicing at these two target genes (Additional file 1: Fig. S2).
PARP1–chromatin structure influences RNA polymerase elongation
We next tested whether PARP1 regulates splicing through regulation of the rate of RNAPII elongation. For this, we used our genome-wide data of PARP1 nucleosome occupancy (GSE56120) in Drosophila S2 cells with PRO-seq data [28] (GSE42117) of transcriptionally engaged RNA polymerase. Analyses of these binding profiles showed that PARP1 and engaged RNAPII are in close proximity within gene bodies. Indeed, we observe a shift ~ 25 bp downstream of the PARP1 signal relative to the RNAPII signal (shown as metagene plots in Additional file 1: Fig. S3). These data suggest that PARP1 may be involved in RNAPII elongation stalling. Next, we investigated if these binding profiles are true in our genes of interest above.
The processivity of RNAPII depends on the phosphorylation state of its carboxy terminal domain (CTD). In particular, the transition between initiation-pausing and productive elongation is marked by phosphorylation on Ser5 and Ser2, respectively [29,30,31]. We therefore asked whether PARP1 influences the recruitment of different forms of phosphorylated RNAPII to exonic regions of our target genes, AKAP and CAPER. Three antibodies that bind specific phosphorylation states of RNAPII were used: (1) 4H8, which recognizes Ser5 phosphorylation (marks initiating and first regions of the gene)—hence referred to as Ser5; (2) H5, which recognizes Ser2 phosphorylation is found mainly in the gene body and toward the 3′ end of the gene. This form is known as the transcriptionally engaged and elongating form of RNAPII—hence referred to as Ser2P; (3) 8WG16, which recognizes hypo-phosphorylated RNAPII found at pre-initiation sites. The 8WG16 antibody has been reported to sometimes also recognize Ser5-phosphorylated RNAPII and could be used to determine total RNAPII. We therefore performed ChIP-qPCR using these antibodies on cross-linked chromatin from wild type (WT) and PARP1 knockdown (KD) cells and analyzed the occupancy of the various forms of RNAPII on the two PARP1 target genes—AKAP and CAPER (Fig. 2 and Additional file 1: Fig: S4). To better assess the correlation between PARP1 reduction and RNAPII occupancy, we calculated the ratio between the occupancies of these polymerase forms and PARP1 at three locations: (1) the immediate preceding constitutive exon; (2) the intervening intron; (3) the proceeding alternative exon of these genes as shown in Fig. 2a. We call this the ‘travelled’ index as it determines the ratio of PARP1 (or RNAPII) occupancy at the proceeding intron or alternative exon, relative to the preceding constitutive exon.
To begin, the occupancy of PARP1 was measured. In WT cells, at the AKAP gene, a 22% reduction in the relative occupancy of PARP1 was observed at the alternative exon 5 compared to constitutive exon 4 (Fig. 2b and Additional file 1: Fig. S4A—blue bar). When using the ‘travelled’ index, this represented a steady decline of PARP1 occupancy from the 5′ constitutive exon 4 to proceeding intron 4–5 and was lowest at the 3′ alternative exon 5 (Fig. 2b). In PARP1 depleted conditions (PARP1 KD cells), a reduction of ~ 80% of PARP1 occupancy was measured at alternative exon 5 relative to constitutive exon 4 (Additional file 1: Fig. S4A—red bar). With the travelled index incorporated, this change represented a further decrease of ~ 45% and ~ 80% in PARP1 occupancy at the proceeding intron 4–5 and at alternative exon 5 relative to the constitutive exon (Fig. 4b—red line). A similar trend was observed at the CAPER gene though the decrease in the relative amount of PARP1 occupancy from constitutive exon 3 to alternative exon 4 was less pronounced. As expected, there was an additional decrease in PARP1 occupancy at all the tested sites in KD cells (Fig. 2c and Additional file 1: Fig. S4A).
Next, we asked whether the observed changes in PARP1 occupancy correlate with changes in the occupancy of RNAPII forms. Using the Ser2P antibody, which recognizes the transcriptionally engaged, elongating form of RNAPII, we showed that depletion of PARP1 correlated with depletion of this transcriptionally engaged RNAPII at the studied exons. These data corroborated our directly measured results comparing alternative versus constitutive exons and also those measured by the travelled index. In WT cells, as measured through direct comparison of occupancy at alternative exon versus constitutive exon, the AKAP gene showed a decrease of ~ 20% (Additional file 1: Fig. S4B). Further supporting these data, the travelled index (Fig. 2d) showed a steady decline of Ser2P from the 5′ constitutive exon toward the 3′ alternative exon. These results were very similar for the CAPER gene (Fig. 2e and Additional file 1: Fig. S4B). Interestingly, in PARP1 knockdown cells, Ser2P decreased even further (~ 40% at intron 4–5 and alternative exon 5 compared to constitutive exon 4 at AKAP) (Fig. 2d and Additional file 1: Fig. S4B), while this decrease in Ser2P was ~ 50% at CAPER gene regions (Fig. 2e and Additional file 1: Fig. S4B). These findings correlated with the observed decrease in PARP1 occupancy at these sites (Fig. 2b, c; Additional file 1: Fig. S4A). Next, we asked if this correlation of PARP1 occupancy is specific for only the elongating Ser2P occupancy. For this, we tested the occupancy of the other phosphorylated forms of RNAPII—the pre-initiating form of RNAPII, also known as hypo-phosphorylated RNAPII (8WG16), and Ser 5 (4H8). In WT cells, the occupancy of 8WG16 at alternative exons over constitutive exons was reduced by 50% for both AKAP and CAPER genes, respectively (Additional file 1: Fig. S4C—blue bars). On the other hand, in KD cells, we observed an increase in the presence of this polymerase form (Additional file 1: Fig. S4C—red bars). Profiling of 4H8, which measures Ser5 which is found at TSSs and gene bodies, showed a slight increase at the alternative exon 5 of AKAP compared to the constitutive exon 4. For CAPER, we observed a large and significant decrease of ~ 80% occupancy of Ser5 (Additional file 1: Fig. S4D—blue bars). In KD cells, AKAP exhibited an increase in this polymerase form while CAPER showed no significant difference in occupancy compared to WT cells (Additional file 1: Fig. S4D—red bars). Finally, analysis of the occupancies of PARP1 and the various RNAPII were recapitulated in cells treated with a second siRNA (siRNA2), thus confirming the effect of PARP1 on the occupancy of RNAPII. In contrast, inhibition of PARylation showed no differences in PARP1 or in the occupancies of the different RNAPII forms when compared to WT cells (Additional file 1: Fig. S5). These data therefore show that PARP1 occupancy and not its PARylation activity exerts an effect on the occupancy of elongating polymerase.
PARP1 depletion disrupts chromatin state and structure
Chromatin context can affect the rate of RNAPII elongation, which in turn, would regulate alternative splicing. After confirming the relationship between PARP1 and RNAPII pausing, we next investigated the type of chromatin context mediated by PARP1 at these alternative exon sites. For this, we mapped the nucleosomes across AKAP and CAPER genes using the nucleosome walking method [32]—a low-resolution technique, which allows gene-specific high-resolution mapping of nucleosome positions along a stretch of DNA [33]. Chromatin is digested with micrococcal nuclease (MNase) to yield mostly mononucleosomal fragments and is then subjected to quantitative real-time PCR (qRT-PCR—see Methods). First, we predicted nucleosome locations based on sequences alone [34] (Fig. 3a, b—top panels). Then, primers were designed to amplify about 80–100-bp-sized amplicons that overlapped by 20–40 bp, tiling across the selected loci of the AKAP and CAPER genes. In this technique, amplification of a product indicates the presence of a protected mononucleosome, while the lack of amplification signifies open chromatin susceptible to MNase digestion. Nucleosome positions and strength of nucleosome occupancy were then calculated using the fold change between MNase-treated samples and undigested genomic DNA at an equivalent DNA input concentration (see Methods).
Based on this method, stable and highly reproducible profiles of nucleosomes were observed across the target gene region (Fig. 3), with clear differences in the nucleosome profile between WT and PARP1 KD cells. For the AKAP gene, two clear observations were made: (1) There is a strong nucleosome in PARP1 KD cells, mapped by primer A8, which was previously absent in WT cells. (2) At the genomic locations mapped by primers A13 to A15, a reduction in nucleosome occupancy in PARP1 KD cells was seen, with a shift of the nucleosomes toward the A16 position (Fig. 3a). For the CAPER gene, we observed a reduction in nucleosome occupancy just before the alternate exon 4 (mapped by primers C8–C10) and in the region containing the alternative exon 4 (mapped by primers C12–C16) (Fig. 3b). Here too, PARP1 knockdown using siRNA2 produced similar results for nucleosome repositioning as seen in cells transfected with PARP1 siRNA1, while PARylation-inhibited cells exhibited no changes in nucleosome positioning when compared to WT cells (Additional file 1: Fig. S6). Since these PARP1-mediated nucleosome rearrangements occur right before the alternate exon, we posit that PARP1 maintains a chromatin structure that would be amenable to transcription elongation by RNAPII in the absence of PARP1.
PARP1 occupancy displays interplay of selective acquisition of histone methylation at genic regions
Several histone marks have been implicated in alternative splicing [20]. Given that our previous data showed interplay between PARP1 and certain histone marks [5], we then sought to determine whether interplay of PARP1 with specific histones could drive the observed chromatin rearrangement and transcriptional elongation machinery at the studied regions. Typically, paused gene regions are marked by bivalent histone marks. In view of this, we used ChIP-qRT-PCR, to measure the occupancy for both the activating mark, H3K4me3, and repressive mark, H3K27me3, in regions that showed the most change in nucleosome structure—regions mapped by primers A14 for AKAP and C8 for CAPER and its surrounding exons (Fig. 3). In WT cells, at both genes, H3K4me3 and H3K27me3 were found at all sites tested (preceding constitutive exon, intervening intron, and alternative exon) with varying levels (Fig. 4). The presence of both the activating H3K4me3 and repressive mark H3K27me3 is indicative of a poised gene region. To further assess whether there is interplay between these histone marks and PARP1, we investigated whether their occupancy is changed in the absence of PARP1. Knockdown of PARP1 resolved this bivalency of H3K4me3 and H3K27me3 marks to H3K4me3. At both of these genes, there was an increase in H3K4me3 occupancy (Fig. 4a for AKAP and 4c for CAPER) and a decrease in H3K27me3 occupancy (Fig. 4b for AKAP and 4D for CAPER). The resultant net gain of H3K4me3—an active histone mark—possibly opens up the chromatin structure allowing for the passage of transcription machinery. These data were recapitulated in cells treated with siRNA2 (Additional file 1: Fig. S7—red bars vs. blue bars), while PARylation-inhibited cells showed no difference relative to WT cells (Additional file 1: Fig. S7—green bars vs. blue bars). In summary, our data are consistent with a model in which binding of the PARP1 mediates or is mediated by specific histone modifications. Additionally, the effect of PARP1 on chromatin structure (nucleosome positioning and occupancy of histone PTMs) is instigated by the direct presence of PARP1 and not its PARylation activity. Thus, at the sites of alternative splicing, PARP1 could play a dual role in stimulating the release of RNAPII pausing and recruiting chromatin modifications that facilitate its release from the paused state.
PARP1 influences RNAPII elongation genome-wide
The observations of the direct effect of PARP1 binding on chromatin structure and histone modification occupancy prompted us to ask whether PARP1 could influence the release of RNAPII from pause sites. We used a modified 3′NT-seq method [35] and the NET-seq method [36] to map the positions of elongating and arrested RNAPII complexes at nucleotide resolution (Fig. 5a) in the presence and absence of PARP1. The 3′NT method effectively isolates mRNA from RNA-RNAPII–chromatin complexes. The presence of the m7G cap on RNAPII transcripts within the first 20–30 nucleotides of transcription allows for the immunoprecipitation of nascent mRNAs using the GFP-elF4E protein (which binds to the m7G cap) bound on magnetic beads. These captured mRNAs were then eluted from beads, purified, and ligated to the 3′ adapter used in Illumina sequencing. Next, the RNAs were fragmented to not only decap the captured 5′ ends of the mRNAs, but also to reduce the size of the mRNAs to ~ 35–100 nucleotides. Fragmented mRNAs were then size selected on a denaturing gel. Next, the 5′ phosphate groups were removed enzymatically, allowing for the ligation of the 5′ sequencing adapter. cDNA libraries were then prepared, and limited PCR amplification with primers that bind to both the 5′ and 3′ adapters was performed. This allowed for the capture and sequencing of only mRNA fragments with both the 5′ end (end after fragmentation) and the original 3′ end of the nascent mRNAs. After PCR, samples were size selected on a 3.3% NuSieve agarose gel. These fragments were treated and analyzed separately, gel purified, and subjected to Illumina high-throughput 50 bp paired-end sequencing. Two biological replicates for WT and PARP1-KD cells were sequenced, generating ~ 30–100 nt reads for each 3′NT-seq sample (LW1 and its corresponding HW1; LW2 and its corresponding HW2). A total of 780 million reads were sequenced, 101 million of which mapped uniquely (i.e., after removing multi-mapped reads and potential PCR duplicates) to the Drosophila genome (Dm6) after additional filtering steps.
In 3′NT-seq experiments, the sequenced read density reflects the abundance of the transcript and the 3′ ends of the nascent mRNAs map the RNAPII position at nucleotide resolution. Thus, assuming there is no degradation, these sequenced RNAs would have the captured 3′ ends of the elongating polymerase just before transcription elongation inhibition by alpha-amanitin. In fact, the resolution afforded by 3′NT-seq and the coverage obtained should provide an in-depth view of genome-wide transcriptional activity. Thus, to begin our analyses, we first compared the reproducibility of the biological replicates using the multiBamSummary tool from DeepTools 2.0 [37]. The biological replicate libraries show strong agreement (Pearson’s coefficient > 0.988), which documents the robustness of our approach. We then compared the sequencing reads normalized by reads per kilobase of transcript, per million mapped reads (RPKM) between WT and PARP1 KD cells. To determine if transcription of specific segments of mRNA genes is targeted by PARP1, we calculated the number of normalized reads for all mRNA genes and divided them into five separate regions: 1000 bases immediately upstream of the start codon (upstream); transcript (from transcription start sites (TSS) to transcription end sites); first 100 bp of the transcript; last 100 bases of the transcript; and 1000 bases downstream of the transcript end (downstream).
A total of 100 bins were created for each region. In the initial global test, using a cutoff of p < 0.01, subtle differences between WT and PARP1 KD samples were observed in the bins for the upstream and transcript regions (Additional file 1: Fig. S8A, B). On the other hand, some differences were observed at downstream regions as well as at the early (first 100 nt) and late gene bodies (last 100 nt) due to PARP1 knockdown (Additional file 1: Fig. S8C–E). Next, using the difference in the percentage of reads occurring in the first 50 bins, we filtered the differentially expressed regions into three groups: Group 1: genes with reads that were shortened; Group 2: genes with reads that were lengthened; Group 3: genes with reads that had differential patterns not specifically related with shortening or lengthening. This analysis resulted in a total of 1786 genes, of which 348 (Group 1) and 307 (Group 2) had lengthened and shortened transcripts, respectively, in KD versus WT cells, and the rest were placed in Group 3. We show the top 20 shortened (Additional file 1: Table S1) and lengthened genes (Additional file 1: Table S2). Examples of shortened and lengthened genes are shown in Additional file 1: Fig. S9A and B, respectively.
To better understand the pattern of these lengthened and shortened genes, we also performed metagene analysis of the 3′NT-seq at different genomic regions. We observed a shifting of RNAPII locations, evidenced by average RNAPII densities in the PARP1 KD, which were substantial, often decreasing (shortened) or increasing (lengthened) (Fig. 5b & c), when compared to the same sites in WT cells. In the KD condition, when transcripts have shortened RNAPII profiles, there is a slight shift in the peak of RNAPII location toward the 5′ end of the gene relative to WT (Fig. 5b). For lengthened genes, one main peak (location) of RNAPII was observed in WT cells. In PARP1 depleted cells, there was a decrease in this peak with a concurrent appearance/increase of a new peak located 3′ of this peak (Fig. 5c), resulting in two prominent peaks.
At the upstream regions, in genes with a shortened RNAPII profile, one peak was detected with little to no change in PARP1 KD cells (Additional file 1: Fig. S10A). As for the lengthened genes, several peaks were present. In cells with depleted PARP1, there was a reduction in the peaks at the 5′ region with a concurrent increase in the peaks at the 3′ end (Additional file 1: Fig. S10B). A similar situation occurred for the RNAPII downstream regions with shortened genes (Additional file 1: Fig. S10C). On the other hand, slight differences were detected with lengthened genes, including a stronger peak emerging toward the 5′ of the WT RNAPII positions (Additional file 1: Fig. S10D). At the other gene regions, the first and last 100 bp gene regions, a shortened RNAPII profile can be seen, similar to the downstream regions (Additional file 1: Fig. S10E and G, respectively). With the lengthened genes, there seems to be a flattening and merging of the two RNAPII peaks in PARP1 KD cells for the RNAPII locations within the first and last 100 bp genic regions (Additional file 1: Fig. S10F and H). According to other studies, a fast polymerase is typically associated with an overall flattening of the RNAPII profile in the termination zone replacing the clear drop-off in RNAPII density that occurs in WT cells. Interestingly, such a flattened profile has been associated with reduced pausing [38]. Finally, we asked whether these genes with shortened or lengthened transcript regions are involved in any functional pathways. For this we performed gene ontology: biological process analysis (GO:BP) using categoryCompare [39]. GO:BP analysis results showed that genes with a shortened RNAPII profile were significantly enriched in genes associated with cell processing (Fig. 5d and Additional file 1: Table S3). On the other hand, the genes with a lengthened RNAPII profile were significantly involved with organismal organization (Fig. 5d and Additional file 1: Table S4).
We also provide NET-seq analyses on our candidate genes, AKAP and CAPER. Prominent peaks of RNAPII density are seen at specific nucleotides or within narrow regions, possibly indicative of pause sites (locations where RNAPII is detected with high probability). For AKAP (Fig. 6a), NET-seq showed a prominent peak at the 5′ of the region, with more RNAPII footprints downstream. Of interest, in KD cells, the main 5′ peak is absent, with more RNAPII footprints downstream compared to WT. Indeed, the last nucleotide added is further along the gene. These results suggest that after the PARP1-mediated block is relieved (through knockdown of PARP1), more RNAPII was able to move downstream to selectively include specific exons over others (Fig. 6a). The situation at the CAPER gene is slightly different. NET-seq in KD cells shows an increase in RNAPII footprint within the preceding intronic region (left panel—Fig. 6b), and a slowing down of RNAPII around the proceeding exonic region located at 7,035,750 bp (Fig. 6b). These results of a fast RNAPII elongation at upstream intronic regions and slow down of RNAPII elongation at the proceeding exons, could explain the selective exon skipping seen in PARP1 KD cells (Fig. 1b). Additionally, these results resonate with previous studies showing that slowing down of RNAPII not only results in exon inclusion, but has been implicated in exon skipping as well [40]. In summary, our analyses for specific genes show a clear RNAPII pausing defect (positive and negative) due to PARP1 depletion. Interestingly, these changes in RNAPII elongation occur at the same location of PARP1-mediated chromatin changes (Figs. 3, 4 and Additional file 1: Figs. S5 and S6), thus supporting our hypothesis that PARP1-mediated chromatin structural rearrangement regulates RNAPII elongation and splicing decisions.