Differential nucleosome occupancy modulates alternative splicing in Arabidopsis

Alternative splicing (AS) is a major gene regulatory mechanism in plants. Recent evidence supports co-transcriptional splicing in plants, hence the chromatin state can impact AS. However, how dynamic changes in the chromatin state such as nucleosome occupancy influence the cold-induced AS remains poorly understood. Here, we generated transcriptome (RNA-Seq) and nucleosome positioning (MNase-Seq) data for Arabidopsis thaliana to understand how nucleosome positioning modulates cold-induced AS. Our results show that characteristic nucleosome occupancy levels are strongly associated with the type and abundance of various AS events under normal and cold temperature conditions in Arabidopsis. Intriguingly, exitrons, alternatively spliced internal regions of protein-coding exons, exhibit distinctive nucleosome positioning pattern compared to other alternatively spliced regions. Likewise, nucleosome patterns differ between exitrons and retained introns pointing to their distinct regulation. Collectively, our data show that characteristic changes in nucleosome positioning modulate AS in plants in response to cold.


Introduction
Plants employ different strategies to control their transcriptional program during the daily cycles of light-dark and in response to environmental stress to confer adaptive responses (Zhu, 2016;Laloum et al., 2017;Lämke & Bäurle, 2017).Recent evidence shows that alternative splicing (AS) regulation is a key gene regulatory mechanism in plants (Calixto et al., 2018;Filichkin et al., 2018;Jabre et al., 2019).In plants and animals, AS is regulated co-transcriptionally (Brody et al., 2011;Tilgner et al., 2012;Li et al., 2020;Zhu et al., 2020).RNA polymerase II (RNAPII) speed during transcription may be affected by the chromatin state that in turn determines AS outcomes (Alexander et al., 2010;Ullah et al., 2018;Zhu et al., 2018).Emerging evidence shows that the chromatin environment has a strong bearing on the splicing process by modulating RNAPII processivity and splicing factors (SFs) recruitment (Nojima et al., 2018;Jabre et al., 2019;Kindgren et al., 2019;Li et al., 2020;Yu et al., 2019;Zhu et al., 2020).Recent native elongating transcript sequencing (NET-Seq) and global run-on sequencing (GRO-Seq) studies from mammals and Arabidopsis thaliana (hereafter Arabidopsis) show that phosphorylation of RNAPII Cterminal domain mediates interactions with the spliceosome and that RNAPII accumulation is associated with different chromatin states (Nojima et al., 2018;Zhu et al., 2018).Remarkably, sequencing of the chromatin-bound nascent RNAs in Arabidopsis revealed that almost all introns are spliced co-transcriptionally and the efficiency of intron removal is more robust in proteincoding genes than noncoding RNAs (Li et al., 2020).Furthermore, it was also demonstrated that co-transcriptional splicing (CTS) efficiency is dependent on the number of exons but not gene length (Zhu et al., 2020).Therefore, appropriate exon-intron definition may be important for CTS in Arabidopsis (Li et al., 2020).For example, nucleosome occupancy was found to be higher on exons and also accompanied by a higher level of RNAPII (Chodavarapu et al., 2010).Therefore, higher nucleosome occupancy on exons is intrinsically associated with exon definition, RNAPII processivity and splicing kinetics (Jabre et al., 2019).Previously, it has been demonstrated that change in RNAPII speed in both directions can influence splicing factor recruitment and splicing efficiency (Dujardin et al., 2014;Godoy Herz et al., 2019;Leng et al., 2020).Therefore, nucleosome positioning may modulate different AS events and their ratios under variable growth and/or stress conditions to alter intron/exon boundaries and provide a context through which AS patterns could be modulated (Jabre et al., 2019).For example, RNAi lines of a chromatin remodeler gene (ZmCHB101) in maize showed altered nucleosome density, RNAPII elongation

Accepted Article
This article is protected by copyright.All rights reserved rate, and changes in splicing patterns under osmotic stress (Yu et al., 2019).Similarly, widespread nucleosome remodelling in rice, as a result of phosphate starvation and cold stress, was associated with differential gene expression (Roy et al., 2014;Zhang et al., 2018).Since cold can influence the RNAPII elongation kinetics in Arabidopsis (Kindgren et al., 2019), we reasoned that rapid cold-induced alternative splicing response in Arabidopsis (Calixto et al., 2018) may be associated with nucleosome remodelling.Henceforth, we used cold treatment as a system of choice to investigate whether it could modulate nucleosome positioning and influence AS.
We performed RNA-Seq and micrococcal nuclease sequencing (MNase-Seq) for Columbia wild type (Col-0) accession of Arabidopsis growing at normal temperature (22 o C) and under cold stress (4 o C) for 24 hours.We observed genome-wide changes in AS and gene expression between Col-0_22 o C and Col-0_4 o C plants.Our results show that temperature-dependent differences in nucleosome positioning are sufficient to modulate different types of AS events and their abundance.Remarkably, exitrons, alternatively spliced internal regions of protein-coding exons (Marquez et al., 2015), can also be distinguished from flanking exons by distinctive nucleosome occupancy to facilitate their recognition by the splicing machinery.

Materials and Methods
The detailed experimental procedure is provided in the supplementary information file (Supporting Information Method S1).Briefly, leaf tissues were harvested from three weeks old Col-0 plants grown at 22°C and cold treated (4°C) for 24 h.Total RNA and nucleosome bound genomic DNA (gDNA) were extracted for Illumina paired-end sequencing using RNA extraction and Qiagen DNA kits, respectively.The raw reads generated from RNA-Seq and MNase-Seq experiments were quality checked using Trimmomatic (Bolger et al., 2014).The high quality reads from RNA-Seq experiment were used to quantify the transcripts expression using Salmon v0.82 (Patro et al., 2017) and AtRTD2-QUASI (Zhang et al., 2017) as reference.Differential Expressed Genes (DEGs) and differential alternatively spliced (DAS) genes were identified using 3D-RNA-Seq pipeline as described previously by (Calixto et al., 2018;Guo et al., 2019).Gene functional enrichment analysis was performed using DAVID v6.8 (Huang et al., 2009a,b).The gene ontology (GO) terms were assigned to DEGs and DAS genes with FDR ≤0.05.AS events, AS event inclusion level (Percent Spliced In -PSI indicates how efficiently sequences of interest are spliced into transcripts) and the difference in this inclusion ( PSI) between Col-0 grown at Δ

Accepted Article
This article is protected by copyright.All rights reserved 22 o C and 4 o C were identified using SUPPA2 (Alamancos et al., 2015;Trincado et al., 2018).Only AS events having a p-value ≤0.05 are identified as DAS events.For MNase-Seq high quality reads were mapped to the TAIR10 Arabidopsis reference genome using Bowtie v1.2.2 (Langmead et al., 2009) with "-m" set to 1 to output only uniquely mapped reads.Improved Nucleosome-Positioning Algorithm (iNPS) was used for accurate genome-wide nucleosome positioning as described previously by (Chen et al., 2014).Differential Nucleosome Positioning (DNP) analysis was performed using DANPOS v2.1.2(Chen et al., 2013).Nucleosome signals were plotted around different genomic regions using deepTools v3.5.0 (Ramírez et al., 2014).

Cold-regulated DE and DAS genes affect different biological processes
In Arabidopsis, nucleosome positioning differentially marks promoter regions, gene bodies as well as exons and introns, indicating a potential link of chromatin architecture to gene expression and splicing regulation (Chodavarapu et al., 2010).To investigate if cold-induced AS in Arabidopsis is regulated by nucleosome occupancy, we performed RNA-Seq and MNase-Seq of Col-0 lines before and after a shift from 22°C to 4°C for 24 hours.Using the previously published 3D-RNA-Seq pipeline (Calixto et al., 2018;Guo et al., 2019) (Supporting Information Method S1), we identified 6252 DEGs and 2283 DAS genes.Of the 6252 DEGs, 3323 were up-regulated and 2929 were down-regulated (Table S1).We observed that most transcriptional changes are associated with genes that do not display splicing changes (70.5% DE only genes).Similarly, a large proportion of splicing changes occur in genes that are not differentially expressed (19.3% DAS only genes).The overlap of DEGs and DAS genes is significant (10.2%; hypergeometric test, p <1.222e -39 ) (Fig. 1a), suggesting that cold stress modulates both transcriptional and AS responses of some genes, which is in line with previously published reports from Arabidopsis (Calixto et al., 2018).Gene functional enrichment analysis of DEGs showed significant (FDR <0.05) enrichment in diverse biological functions including circadian rhythm, cold stress, and photosynthesis regulation.Cellular components terms enrichment for DEGs were mainly for plasma membrane, and vacuole (Fig. S1a).DAS genes showed significant (FDR <0.05) enrichment in mRNA processing, RNA splicing and protein phosphorylation, whereas those of cellular components and molecular functions were mainly enriched in the nucleus and mRNA/ATP/protein binding

Nucleosome occupancy modulates variety and abundance of alternative splicing events
To investigate how nucleosome occupancy modulates splicing, we performed nucleosome positioning and DNP analysis (Supporting Information Method S1).We detected 19233 significant Differential Positioned Nucleosomes (DPNs) upon shifting plants from 22 o C to 4 o C for 24 hours that were significantly (Supporting Information Method S1) associated with 7357 genes (Tables S4, S5, and S6).Interestingly, 833 (9.5%) (hypergeometric test, p < 1.498e -24 ) coldinduced DAS genes displayed changes in nucleosome occupancy (Fig. 1c).We first profiled nucleosome occupancy across exons and flanking regions, where we could detect a significant drop of nucleosome occupancy signals at 4°C around exons and flanking regions (one-tailed t test, p <0.0001) (Fig. S3).Nucleosome occupancy is strongly associated with negative or positive AS regulation Next, we interrogated how nucleosome occupancy levels, for the same set of genes, differ between DAS and non-DAS genes under normal and cold conditions.For that, we profiled nucleosome occupancy levels across the exons of DAS and non-DAS genes.We found relatively lower nucleosome occupancy for DAS compared to non-DAS exons at 22 o C and as well as 4 o C (Onetailed t test, p <0.0001, Fig. 3a).Since nucleosome occupancy globally drops under cold conditions, we sought to investigate how nucleosome occupancy would correlate with splice junctions affected differently by cold stress.Therefore, we grouped the SJs of the AS events (pvalue < 0.05) obtained from SUPPA based on their ΔPSI value to obtain positively, negatively, and unaffected SJs (Supporting Information Method S1).Interestingly, cold stress positively regulates 1208 SJs, negatively regulates 1054, and leaves 673 SJs un-affected (Fig. 3b; Table S7).
This data strongly supports previous data showing that cold-stress induced AS in plants (Calixto et al., 2018).To profile nucleosome occupancy around these SJs, we plotted nucleosome density of negatively-, positively-affected, and unaffected SJs for Col-0 grown at 22 o C and 4 o C (Fig. 3c).
Remarkably, nucleosome profiles of unaffected SJs for both, the donor and acceptor site are significantly different compared to negatively or positively affected SJs at both temperatures (One-way ANOVA, p <0.01, Fig. 3C).Additionally, we also detected a significant association between negatively and positively affected SJs with regions associated with DNPs (Fisher exact

Accepted Article
This article is protected by copyright.All rights reserved test, p-value <0.001).Overall, our results show that changes in nucleosome occupancy levels across intron-exon junctions and exons are likely to regulate splice site selection and subsequently modulate splicing regulation in both positive and negative manner.

Characteristic nucleosome occupancy patterns define exitrons
EIs have a lower guanine-cytosine (GC) content than adjacent sequences of EI-containing exons (Marquez et al., 2015).Therefore, we asked whether differential GC content in EI sequences is associated with nucleosome occupancy to distinguish them from flanking exonic regions.To answer this, we profiled nucleosome occupancy across ~2400 EIs identified in Arabidopsis (Marquez et al., 2015;Zhang et al., 2017) and 500 bp upstream and downstream from their starts and ends.We found sharp peaks of nucleosome occupancy located before the start and after the end of EIs and slightly lower occupancy in the middle of EIs, which is different from nucleosome patterns observed across exons (Fig. S3).Additionally, we detected a decrease of nucleosome occupancy across exitrons under cold stress (One-tailed t test, p <0.0001, Fig. 4a).Comparison of nucleosome occupancy levels over EIs grouped into four bins according to the PSI values showed that higher exitron inclusion correlates with higher nucleosome occupancy and showed variable levels under normal and cold conditions (Fig. 4b), hence pointing towards their regulation under cold stress.This pattern differs from the one observed for IRs, where IRs with higher inclusion levels display less nucleosome occupancy (Fig. 2b).Interestingly, in this respect, EIs are more similar to ES, A5'SS and A3'SS events due to their exonic features (Fig. 2b).Since EIs have a higher GC content than IRs and constitutive introns (Marquez et al., 2015), we compared their nucleosome profiles in normal and cold conditions.We observed that nucleosome occupancy levels are higher for EIs compared to IRs at 22 o C and 4 o C (One-way ANOVA, p <0.0001, Fig. 4c).This implies that chromatin structure plays different roles in the definition and splicing of IRs and EIs.Overall, these results revealed the importance of nucleosome occupancy in defining EIs and their distinction from IRs to regulate their AS profiles under normal and cold conditions.

Discussion
Recent evidence from Arabidopsis shows that transcription and splicing are largely coupled (Dolata et al., 2015;Hetzel et al., 2016;Ullah et al., 2018;Jabre et al., 2019;Jia et al., 2020), and that epigenetic features in plants regulate transcriptional activity and differentially mark exons and

Accepted Article
This article is protected by copyright.All rights reserved introns (Zhu et al., 2018).Not surprisingly, very recent studies employing sequencing of chromatin-bound RNAs reveal that almost all introns in Arabidopsis are spliced cotranscriptionally (Li et al., 2020).Furthermore, RNAPII elongation speed has been found to be slower in nucleosome-rich exons allowing more time for the splicing process to take place (Chodavarapu et al., 2010;Zhu et al., 2018).However, how the chromatin environment influences different types of AS events and their ratios under variable growth and stress conditions remains elusive in plants.Since splicing/AS regulation is achieved by the context of the cis-regulatory sequences as well as the chromatin environment (Reddy et al., 2013), it is important to understand the relative contributions of epigenetic landscapes.In this study, using Arabidopsis Col-0 ecotype  (Iannone et al., 2015).Since histone modification modulate AS in humans (Luco et al., 2010), similar mechanism also may modulate splicing variation in plants.Intriguingly, exitrons also show distinctive nucleosome occupancy (Fig. 4), which may help to differentiate them from flanking exonic regions.Furthermore, despite their classification as a group of IR, exitrons display different nucleosome patterns compared to retained introns; pointing to their distinct regulation.
We propose that these changes in nucleosome occupancy may provide the basic definition to exons and introns to coordinate RNAPII processivity.However, it is apparent that the splicing process is also fine-tuned by various trans-regulatory factors and histone modifications under variable growth and stress conditions (Kindgren et al., 2019;Zhu et al., 2020).Our data support this notion, and it is likely that higher nucleosome occupancy may regulate RNAPII accumulation around splice sites and enable SFs recruitment to facilitate and/or modulate splicing variation.
Interestingly, RNAPII elongation speed in Arabidopsis would be much slower after clearing a

Accepted Article
This article is protected by copyright.All rights reserved 3'SS and towards the end of an exon, and may not provide sufficient time (because of higher speed in plant introns) for RNAPII to recognise the 5'SS (Kindgren et al., 2019).Furthermore, recent findings show that RNAPII accumulates upstream of the 5'SS, potentially to provide additional time/checkpoint to regulate splicing in mammals and plants (Kindgren et al., 2019;Nojima et al., 2015).Therefore, variation in nucleosome occupancy with an additional peak just after 5'SS (Fig. 3a) may mediate RNAPII accumulation and influence CTS (Nojima et al., 2015;Kindgren et al., 2019).Arguably, this is why 5'SS splicing dynamics are much more complicated and the scanning splicing machinery has to travel to the branch point/ polypyrimidine tract to complete lariat formation and process 5'SS.Beggs and colleagues proposed, the initial propensity of splicing is low but increases subsequently to allow accumulation of splicing precursors to improve splicing efficiency in subsequent and/or successive reactions (Aitken et al., 2011).These findings are in broad agreement with CTS in Arabidopsis as the CTS process is more efficient in genes with multiple introns/exons and is independent of the gene length (Zhu et al., 2020).
Mutations at the 3'SS and 5'SS impact transcription initiation and a mutant 3'SS reduces the first step of CTS in yeast (Aitken et al., 2011).Similarly, splicing dynamics of the human beta-globin gene which fails to form lariat formation and complete 5'SS when a deletion removes the polypyrimidine tract and AG dinucleotide at the 3'SS (Reed & Maniatis, 1985).Therefore, it is tempting to speculate that nucleosome occupancy and/or histone decorations may be more important in the 5' regions of exons providing a checkpoint to the elongating RNAPII to help recognise 5'SS, form lariat and cleave at the 5'SS and 3'SS.It is evident that efficient splicing/AS is dependent on an optimum RNAPII elongation speed and any variation (slow or fast) results in changes in splicing patterns in humans and plants (Dujardin et al., 2014;Godoy Herz et al., 2019;Leng et al., 2020).
Collectively our data points towards the importance of epigenetic features such as nucleosome occupancy for plants grown under different and recurrent growth and stress conditions.

Fig. 1
Fig. 1 Cold-induced changes in gene expression, alternative splicing, and nucleosome occupancy Arabidopsis thaliana.(a) Venn diagram displaying the overlap between DEGs and DAS genes.(b) Histogram representing the number of DAS events detected using RNA-seq data upon cold stress.(c) Venn diagram displaying the overlap between DAS genes and the genes detected within the differentially positioned nucleosome regions.DEGs and DAS are differentially expressed and alternatively spliced, genes, respectively.DNPs and DNPs-Genes are differentially positioned nucleosomes and the genes associated with them, respectively.The p-value (hypergeometric test) relates to the significance of overlap.A5'SS: Alternative 5' splice site, A3'SS: Alternative 3' splice site, IR: Intron retention events without exitrons, MX: Mutually exclusive exons, ES: Skipped exon, AF: Alternative first exon, EI: Exitrons, and AL: Alternative last exon.

Fig. 2
Fig. 2 Association of nucleosome occupancy with different AS events and their different ratios in Arabidopsis thaliana.(a) The association of nucleosome occupancy with different AS events.The x-axis is the position relative to the acceptor site (left) and donor site (right); the y-axis is the average of nucleosome signal for the selected genomic regions.ANOVA-test has been performed to detect the significance of differential nucleosome occupancy around acceptor site (p = 0.015) and donor site (p = 0.039) of different AS events of at 22 o C, and the donor (p = 0.0138) and the acceptor sites (p = 0.0196) of different AS events at 4 o C (b) Nucleosome profiles for different types of AS events grouped based on their PSI.ANOVA-test has been performed to detect significance of differential nucleosome occupancy of different PSI groups for at 22 o C and 4 o C, respectively; around the acceptor site of A3'SS (p = 0.0129, p = 0.00112), A5'SS (p = 0.00033, p = 0.0112 ), ES (p = 0.000129, p = 0.00234), and IR (p = 0.00236, p = 0.132), events.The x-axis is the position relative to the acceptor site; the y-axis is the average of nucleosome signal.ES: Exon skipping, A3'SS: Alternative 3'SS, A5'SS: Alternative 5'SS, and IR: Intron retention.Constitutive exons or introns are coloured in yellow, whereas exons/introns involved in the splicing event are coloured in blue.Curved lines indicate a splicing event.Red arrow pointing towards differences in scaling used to plot nucleosome profiles for 22 o C and 4 o C. Blue arrows indicate regions with significant changes in nucleosome occupancy.

Fig. 3
Fig. 3 Profiles of nucleosome occupancy across DAS, non-DAS exons and alternatively spliced junctions in Arabidopsis thaliana.(a) Nucleosome profiles are plotted against the DAS and non- Then, we sought if changes in nucleosome occupancy around the splice site can modulate different AS events.For that, we profiled nucleosome signals of 3032 cold-regulated DAS events that showed significant (p-value <0.05) ΔPSI values upon cold stress.Interestingly, different DAS events displayed significant changes in nucleosome occupancy signals around the donor and acceptor sites of different AS events at 22 o C and 4 o C (One-way ANOVA-test, p <0.01, Fig2a).We also could detect an overall drop of nucleosome occupancy level for all AS events upon cold stress.For example, nucleosome occupancy is relatively higher for ES event at 22 o C compared to 4 o C (One-tailed t test, p <0.0001, Fig2a) potentially impacting exon definition and loss of exons from different transcripts (Fig.2a).Since nucleosome occupancy levels differentially associate with various types of AS events between plants grown at different temperatures, we explore if this relationship holds true to explain the ratios of these AS events.For that, we grouped PSI values (Supporting Information Method S1) for different AS events detected in plants grown at both temperature conditions into four bins and aligned nucleosome peaks 200 base pairs (bp) upstream and downstream the exon (or intron) for which the PSI value was calculated.It is notable that for ES, A3'SS and A5'SS, alternative regions with higher inclusion levels (higher PSI values) display more nucleosome occupancy across the splice sites, whereas IRs with higher inclusion levels display less nucleosome occupancy across the splice sites (One-way ANOVA-test, p <0.001, Fig 2.b).Collectively, the differences in nucleosome occupancy levels detected for plants grown at different temperatures for different AS events or for the same PSI group within the same AS event show that alternative regions involved in different types of AS events may be associated with specific epigenetic features potentially influencing local splicing events and abundance of transcripts.
This article is protected by copyright.All rights reserved sought to plants, we demonstrate that cold-induced DAS is accompanied by changes in nucleosome occupancy levels.Although nucleosome occupancy falls globally under cold conditions in Arabidopsis; nonetheless, nucleosome profiles around intron/exon boundaries among different PSI groups, negatively and positively affected AS events displayed characteristic patterns.Further work is needed to understand how variable nucleosome occupancy modulates RNAPII (Chodavarapu et al., 2010)ve splicing in plants.However, since nucleosome occupancy and RNAPII density has a close relationship in Arabidopsis(Chodavarapu et al., 2010), it is likely that chromatin architecture plays a similar role in plants and animals as progesterone treated breast cancer cells displayed weaker nucleosome densities and lower RNAPII accumulation, resulting in alteration in splice site recognition and exon skipping