The histone variant H2A.Z promotes splicing of weak introns.

Multiple lines of evidence implicate chromatin in the regulation of premessenger RNA (pre-mRNA) splicing. However, the influence of chromatin factors on cotranscriptional splice site usage remains unclear. Here we investigated the function of the highly conserved histone variant H2A.Z in pre-mRNA splicing using the intron-rich model yeast Schizosaccharomyces pombe Using epistatic miniarray profiles (EMAPs) to survey the genetic interaction landscape of the Swr1 nucleosome remodeling complex, which deposits H2A.Z, we uncovered evidence for functional interactions with components of the spliceosome. In support of these genetic connections, splicing-specific microarrays show that H2A.Z and the Swr1 ATPase are required during temperature stress for the efficient splicing of a subset of introns. Notably, affected introns are enriched for H2A.Z occupancy and more likely to contain nonconsensus splice sites. To test the significance of the latter correlation, we mutated the splice sites in an affected intron to consensus and found that this suppressed the requirement for H2A.Z in splicing of that intron. These data suggest that H2A.Z occupancy promotes cotranscriptional splicing of suboptimal introns that may otherwise be discarded via proofreading ATPases. Consistent with this model, we show that overexpression of splicing ATPase Prp16 suppresses both the growth and splicing defects seen in the absence of H2A.Z.

Functional cross-talk between chromatin, transcription, and messenger RNA (mRNA) processing is well documented (de Almeida and Carmo-Fonseca 2014). Multiple lines of evidence support a model in which pre-mRNA splicing is primarily cotranscriptional, and splicing factors are recruited to, and act upon, nascent RNA while it is engaged with the RNA polymerase II (RNAP II) complex (Merkhofer et al. 2014). Chromatin structure is an important foundation for the coordination of splicing and transcription, providing both a recruitment platform for factors in both pathways and an environment to modulate the kinetics of transcription, which influences splice site recognition and the kinetics of splicing (Gomez Acuna et al. 2013). As such, ablating chromatin factors or altering their distribution leads to defects in splicing efficiency and changes in alternative splicing patterns (Naftelberg et al. 2015). However, the mechanisms by which chromatin influences splice site choice are still poorly understood.

Splicing, the process of removing intronic sequences from pre-mRNA, is carried out by the spliceosome, a highly conserved macromolecular machine, in an intricate, stepwise manner (Supplemental Fig. S1; Wahl et al. 2009). Spliceosome assembly and catalysis are tightly regulated to produce both canonical and alternative transcripts from a diverse set of cis-splicing sequences (Ast 2004). The splicing reaction is driven forward by ATPases that control the rate of spliceosomal rearrangements while providing an opportunity for proofreading of the substrates (Staley and Guthrie 1998; Koodathingal and Staley 2013). Proofreading allows the spliceosome to identify introns with suboptimal cis-splicing sequences that may take more time to process or discard incorrect substrates (Staley and Guthrie 1998; Koodathingal and Staley 2013). In organisms with many diverse intronic signals, auxiliary splicing factors aid splice site selection, but precisely how additional gene expression factors, such as chromatin, participate remains unclear (Chen and Manley 2009).

The distribution and modification pattern of nucleosomes not only serves as a marker of gene boundaries but also provides information across the gene body during gene expression (Campos and Reinberg 2009). Nucleosomes, by virtue of their enriched occupancy in exons, are thought to play a role in defining exon–intron boundaries, and differential histone modifications have been shown to influence splicing efficiency and alternative splicing outcomes (Schwartz and Ast 2010). Reciprocally, splicing factors and intron/exon architecture have been shown to direct histone modifications and positioning (de Almeida et al. 2011; Kim et al. 2011; Bieberstein et al. 2012). In addition to histone modification, replacement of canonical histone proteins with histone variants is another way to add specialization to the nucleosome. H2A variant H2A.Z is conserved across 90% of eukaryotes with 60% sequence similarity (Zlatanova and Thakar 2008). Deposition of the H2A.Z/H2B dimer in exchange for a canonical H2A/H2B dimer is mediated by the Swr1 complex (Krogan et al. 2003; Kobor et al. 2004; Mizuguchi et al. 2004) and typically occurs at the +1 and −1 nucleosomes around the transcription start site (TSS) (Subramanian et al. 2015). Recently, H2A.Z has also been found to incorporate into gene bodies in metazoans (Weber et al. 2010; Tolstorukov et al. 2012). Notably, closely related variant histones in Drosophila melanogaster (H2A.V) and Homo sapiens (H2A.Bbd) are specifically enriched at intron–exon boundaries (Weber et al. 2010; Tolstorukov et al. 2012).

Our initial efforts to examine how the chromatin landscape influences splicing used Saccharomyces cerevisiae (Wilmes et al. 2008). Although genome-wide genetic interaction mapping uncovered many connections between splicing and other gene expression pathways, the scope of this analysis was limited by the fact that the budding yeast genome contains only ∼250 introns, and most adhere to strong consensus cis-splicing signals. In order to extend this analysis to a more metazoan-like system, we turned to the fission yeast Schizosaccharomyces pombe, the genome of which encodes >5000 introns that use a diverse array of cis-splicing signals and a spliceosome that shares more conservation with metazoans (Wood et al. 2002; Kuhn and Kaufer 2003; Sridharan et al. 2011). Furthermore, a subset of fission yeast introns undergoes alternative splicing in the form of intron retention in response to certain environmental conditions (Awan et al. 2013; Bitton et al. 2015; Stepankiw et al. 2015). Thus, using the fission yeast system, we can better interrogate how a diversity of intron strengths and splicing decisions may be impacted by chromatin.

Building on work from Patrick et al. (2015), in which a genome-wide genetic screen in S. pombe revealed extensive connections between splicing and chromatin factors, here we identified a role for the Swr1 chromatin remodeling complex and the histone variant H2A.Z in splicing weak introns in S. pombe. Specifically, introns with suboptimal architecture and nonconsensus cis-splicing sequences require H2A.Z for their efficient removal. The magnitude of this effect is enhanced further during temperature stress. Our data show that H2A.Z is enriched in intron-containing genes, located primarily within the +1 nucleosome, near the intron/exon boundaries of weak introns. Overexpression of Prp16, a splicing ATPase with proofreading activity, suppresses the temperature-dependent growth defect and nonconsensus intron retention defect of strains lacking H2A.Z or the Swr1 ATPase. Additionally, mutating nonconsensus splice sites of a model intron to the consensus sequence was sufficient to relieve the requirement for H2A.Z in promoting splicing. Taken together, these data suggest a crucial role for H2A.Z in marking weak introns, allowing for additional inspection by the spliceosome and prevention of spurious discard.

Results

The Swr1 complex has strong negative genetic interactions with the spliceosome

An epistatic miniarray profile (EMAP), a genome-wide genetic screen for interactions of gene products, in S. pombe previously uncovered a number of strong genetic interactions between splicing factors and the Swr1 nucleosome remodeling complex (Patrick et al. 2015). To further understand this observation, we performed a permutation analysis on the S. pombe EMAP to ask which gene ontology (GO) processes and complexes are enriched for genetic interactions with the Swr1 complex. Indeed, we found that “mRNA metabolic process” is among the top hits for this GO process analysis (Fig. 1A). To extend the resolution of the GO complex term “spliceosome,” we manually curated subcomplexes of the spliceosome and included them with the GO term complex analysis (Supplemental Table S1; Wahl et al. 2009; Cvitkovic and Jurica 2013). This permutation analysis of GO complexes identified an enrichment in negative interactions between the Swr1 complex and several early splicing complexes, such as splicing complex A, SF3A, splicing complex B, U1 snRNP, and U2 snRNP. These subcomplex definitions are not mutually exclusive, as SF3A is associated with the U2 snRNP, and both U1 snRNP and U2 snRNP are part of splicing complexes A and B. While we observed interactions across the splicing cycle, the strongest enrichment is with the early assembling spliceosome factors of the U1 snRNP and U2 snRNP, which promote the recognition of the 5′ splice site (5′SS) and branch point (BP) (Supplemental Fig. S1).

Figure 1.

The Swr1 complex has strong genetic interactions across the spliceosome. (A) Processes and complexes enriched for significant negative genetic interactions with the Swr1 complex in the S. pombe EMAP. Processes are defined by the S. pombe GO Slim database (PomBase). Complexes are from the S. pombe cellular component GO list (PomBase), with splicing subcomplexes manually defined as in Supplemental Table S1. Complexes involved in chromatin remodeling are not shown. Significance is defined as Bonferonni corrected P-value < 0.00005 for processes and P-value < 0.005 for complexes. (B) Heat map of genetic interaction S scores between Swr1 complex factors and splicing factors. Splicing factors shown have at least one S score either >2.0 or less than −2.5. Gray squares indicate no data. S. cerevisiae gene names in are paretheses. (C) Serial dilutions (1:5) of wild type (WT), the pht1Δ mutant, splicing factor mutation alleles representing different splicing subcomplexes, and their double mutants with pht1Δ, grown at 16°C, 30°C, and 37°C.

The majority of splicing factors in fission yeast are essential and are perturbed in the EMAP screen as “decreased abundance by mRNA perturbation” (DAmP) alleles, where an antibiotic resistance cassette has been inserted into the gene between the stop codon and the 3′ untranslated region (UTR) to reduce the gene's expression (Schuldiner et al. 2005). However, subsequent reports show that DAmP alleles of RNA-binding proteins can be overexpressed (Lund et al. 2008; Patrick et al. 2015). Given the varying phenotypic strength reported for DAmP alleles, we were encouraged to see that the magnitudes of genetic interaction scores are consistent across double mutants between the Swr1 complex and the splicing factor DAmP strains throughout the spliceosome (Fig. 1B; Supplemental Fig. S2A). Additionally, not all splicing factor DAmP alleles yielded strong genetic interactions (Supplemental Fig. S2B).

To more closely examine the interactions of the Swr1 complex with splicing factors at multiple steps of the splicing cycle, we manually generated double-mutant strains through crosses and performed serial dilution growth assays (Fig. 1C). Focusing on pht1, the gene encoding the variant histone H2A.Z and primary substrate of the Swr1 complex, and swr1, the gene encoding the ATPase and organizational platform of the Swr1 complex, we crossed deletions of each of these chromatin factors to alleles of splicing factors that have roles throughout the splicing cycle, including usp107, a U1 snRNP factor; sap14Δ, a U2 snRNP factor; dim1, a tri-snRNP factor; cwf15, a nineteen complex (NTC) factor; prp22, a DEAH-box ATP required for mRNA release after the second catalytic step of the splicing reaction; and prp43, a DEAH-box ATPase required for spliceosome disassembly after proofreading-induced discard or the completion of the splicing reaction (Ohi et al. 2002; Carnahan et al. 2005; Newo et al. 2007; Koodathingal et al. 2010; Shao et al. 2012; Cvitkovic and Jurica 2013). As predicted by the negative genetic interaction scores in the EMAP, combining the pht1Δ strain with these six splicing factor alleles results in varying levels of synthetic sickness in the double-mutant strains at 30°C. The serial dilution profiles of swr1Δ and splicing factor double-mutant strains were similarly synthetic sick at 30°C (Supplemental Fig. S2C).

The Swr1 complex is required for pre-mRNA splicing of a subset of introns

Given the strong genetic interactions between the spliceosome and the Swr1 complex, we wanted to investigate what molecular changes were responsible for the observed phenotypes. To test whether changes to splicing were contributing to the Swr1 complex/spliceosome double-mutant strain growth defects, we performed genome-wide splicing-specific microarray analysis (Fig. 2A). These microarrays contain oligonucleotide probes that are specific for exon, intron, and spliced exon–exon junction regions for the vast majority of intron-containing genes in the S. pombe genome, allowing us to compare total mRNA, pre-mRNA, and mature mRNA between two strains (Lipp et al. 2015; Patrick et al. 2015). From these individual measurements, we generated a “splicing index” score, which reports the ratio of pre-mRNA to mature mRNA multiplied by the level of total mRNA to normalize changes in transcription. We compared total RNA from an isogenic wild-type strain with a variety of single- and double-mutant strains of the Swr1 complex and the spliceosome.

Figure 2.

H2A.Z is required for pre-mRNA splicing. (A) Cartoon depicting splicing-specific microarray probes and splicing index. (B) Splicing profile of pht1Δ and usp107 single-mutant and pht1Δ usp107 double-mutant strains. The EMAP S score for pht1Δ usp107 was 0.455335 (neutral). cDNA from each single- and double-mutant strain was competitively hybridized on the splicing-specific microarray against that from wild type. A splicing index value was calculated for each intron by normalizing the log2 ratio of the intron change to junction change, multiplied by the exon change. The heat map shows the splicing index score for most introns in S. pombe of the indicated strain compared with wild type. Gene order along the Y-axis was the same for all arrays. The stacked bar plots below show the number of introns with a splicing index score >0.3 or less than −0.3. The color scale shows the distribution of the severity of the splicing defects. (C) Stacked bar plots of splicing index scores for the pht1Δ strain. The method was the same as in B except that cultures were grown to mid-log phase at 30°C and then shifted for 9 h to 16°C where indicated. Data for 30°C are the same as in B. (D) RT-qPCR validation of the splicing defect observed in pht1Δ at 16°C microarray on three introns. Bars represent the ratio of intron to exon signal. Error bars represent the SEM of three biological replicates. (*) P < 0.05; (**) P < 0.01, calculated by one-tailed t-test

In order to understand whether the synthetic sickness that we observed correlated with a molecular splicing defect, we first analyzed the splicing efficiencies of pht1Δ, usp107, and the double-mutant strain pht1Δ usp107 (Fig. 2B). We were particularly interested in a U1 snRNP factor, such as usp107, because of the enrichment of splicing complex A and U1 snRNP in our permutation analysis of complexes (Fig. 1A). usp107 is overexpressed when assayed by RT-qPCR, possibly creating a dominant-negative allele (Supplemental Fig. S3A). The double-mutant pht1Δ usp107 strain had an exacerbated splicing defect beyond either single-mutant strain, with 567 introns retained, suggesting that fission yeast are especially reliant on Pht1 for splicing in the context of an already defective spliceosome. However, because the double-mutant pht1Δ usp107 strain is sicker than the single mutants, we cannot rule out the possibility that poor growth indirectly accounts for the exacerbation of intron accumulation.

These observations led us to speculate that Pht1 and/or Swr1 may be particularly important for splicing when cells are grown under nonideal conditions. Because many strains harboring mutations in splicing factors are sensitive to changes in temperature, likely due to defects in the kinetics of spliceosomal rearrangements, we repeated the above serial dilution growth assays at 16°C and 37°C (Fig. 1C; Rosbash et al. 1981; Noble and Guthrie 1996). Intriguingly, the pht1Δ strain has been described previously as cold-sensitive at 16°C, while the swr1Δ strain has been shown to be heat-sensitive at 37°C (Supplemental Fig. S3B; Ahmed et al. 2007). This observation indicates that H2A.Z may be important for recruiting or stabilizing gene expression complexes such as the spliceosome during thermal stress. Indeed, at 16°C, the double mutants of usp107, sap14Δ, dim1, or prp43 with pht1Δ are synthetic lethal, and cwf15 and prp22 are strongly synthetic sick. The serial dilution profiles of swr1Δ and splicing factor double-mutant strains were similarly exacerbated compared with either single mutant but under heat stress at 37°C (Supplemental Fig. S2C). The observation that growing these double mutants under temperature stress exacerbates growth defects supports the hypothesis that the interplay between the Swr1 complex and the spliceosome is particularly important for splicing when yeast are grown in nonideal kinetic conditions.

To test this hypothesis, we grew the pht1Δ and swr1Δ single-mutant strains under temperature stress of 16°C and 37°C, respectively, and looked for an intron retention defect beyond the level observed during growth at 30°C (Fig. 2C). Indeed, we observed a significant increase in the number of affected introns when pht1Δ cells were grown at 16°C, with 883 introns affected (compared with 37 introns at 30°C). A similar trend was observed for swr1Δ at 37°C, with 155 introns affected compared with four introns at 30°C (Supplemental Fig. S3C). While the loss of the Swr1 ATPase has been shown to dramatically decrease the level of H2A.Z in chromatin, Swr1-independent incorporation and a residual fraction of deposited H2A.Z may account for this moderate intron retention profile in swr1Δ compared with pht1Δ (Papamichos-Chronakis et al. 2011; Sadeghi et al. 2011). We confirmed the intron retention defect in the pht1Δ strain at 16°C by measuring spliced and unspliced RNA from several affected introns using RT-qPCR (Fig. 2D).

To further examine the role of the Swr1 complex in pre-mRNA splicing, we analyzed a published RNA sequencing (RNA-seq) data set of swc5Δ, a deletion strain of a component of the Swr1 complex responsible for the histone replacement activity along with swr1 (Wu et al. 2005; Clement-Ziza et al. 2014). We compared the swc5Δ data sets with wild-type data sets and observed a statistically significant increase in intronic reads normalized to total reads for each data set in the swc5Δ strain, indicating a global increase in retained introns and further supporting a role for the Swr1 complex in promoting pre-mRNA splicing in fission yeast (Supplemental Fig. S3D).

To address the concern that changes in global transcription may be responsible for the intron accumulation defects observed, we used published pht1Δ and swr1Δ expression microarrays and swc5Δ RNA-seq (all three performed at 30°C) to look for splicing factor genes with greater than twofold expression level changes. We did not find consistent overlap of gene expression changes between these three data sets. Additionally, we selected eight splicing factor genes from these published data sets and performed RT-qPCR comparing wild-type with pht1Δ strains at 16°C (Supplemental Fig. S4). We observed no changes in expression in these splicing factor genes greater than twofold in the pht1Δ strain background. Taken together, these data strongly suggest that changes in global transcription are not responsible for the intron accumulation defects that we observed (Kim et al. 2009; Clement-Ziza et al. 2014; Larson et al. 2016).

H2A.Z-affected introns are more likely to be 5′ in the gene and contain nonconsensus cis-splicing sequences

Because we observed intron specificity in our splicing-specific microarray of pht1Δ at 16°C, we asked whether the subset of pht1Δ-affected introns contained cis-splicing signals that made them more sensitive to loss of H2A.Z. To address this question, we analyzed different intron features using a logistic regression model to see whether there was a correlation between cis-splicing signals and intron retention in the pht1Δ strain at 16°C (Lipp et al. 2015). Because H2A.Z is known to be enriched in the promoter regions of genes, we first tested whether the distance of an intron from the TSS correlated with intron retention in the absence of pht1 at 16°C (Fig. 3A; Subramanian et al. 2015). Interestingly, we observed that the likelihood of intron retention in pht1Δ at 16°C was inversely proportional to proximity to the TSS, meaning that introns located in the 3′ end of a gene were less likely to accumulate in pht1Δ cells at 16°C. The increased sensitivity of introns at the 5′ ends of genes to the ablation of pht1 is consistent with the enrichment of H2A.Z in the promoter region of genes.

Figure 3.

Loss of H2A.Z impairs splicing of weak introns at the 5′ ends of genes. (A) Schematic representation of an S. pombe intron with arrows indicating features of intron architecture. The bar graph shows the probability of the indicated feature increasing the likelihood of intron retention, represented as log odds or log odds per nucleotide (positive log odds at a given feature indicate that increased length or frequency correlated with intron retention.) Affected introns were defined as having a corrected P-value of <0.05. Error bars represent the SEM of a logistic regression model (Lipp et al. 2015). (*) P < 0.05; (**) P < 0.01; (***) P < 0.001. (B) Schematic representation of the S. pombe 5′SS, BP, and 3′SS, with consensus nucleotides indicated in black and nonconsensus nucleotides indicated in red. The probability of a specific nucleotide change increasing the likelihood of retention is represented as log odds (positive log odds at a given position indicate that nonconsensus nucleotides correlate with intron retention). Statistics are the same as above. (C) Histogram of 5′SS and BP strength in S. pombe introns, plotted by their distance in base pairs from the TSS. Introns in the bottom decile of 5′SS and BP strength are enriched in the 5′ ends of genes. P-values are from a Kolmogorov-Smirnov (K-S) test.

We next examined correlations between intron retention and other architectural features of introns, including intron length, distance between the BP and the 3′ splice site (3′SS), and presence of a polypyrimidine tract. Previous work in S. pombe has shown that deviation from average length or absence of a polypyrimidine tract is associated with lower splicing efficiency (Sridharan et al. 2011). Indeed, we found that introns are more likely to accumulate in the absence of H2A.Z if they are longer, have a longer distance between the BP and 3′SS, and lack a polypyrimidine tract (Fig. 3A), suggesting that Pht1 is particularly important for splicing pre-mRNAs with suboptimal intron architectures.

To determine whether pht1Δ-affected introns were more or less likely to be enriched for particular cis-splicing sequences, we looked at the 5′SS, BP, and 3′SS of introns that accumulated in pht1Δ cells at 16°C. We based this analysis on the fact that, although cis-splicing sequences in S. pombe are less constrained than those in S. cerevisiae, there is still an overall consensus sequence for each feature (Kuhn and Kaufer 2003; Stepankiw et al. 2015). Our analysis revealed that particular nucleotide deviations from the consensus cis-splicing sequences (5′SS, BP, and 3′SS) are associated with significantly increased probability that an intron will be retained in pht1Δ at 16°C (Fig. 3B). As these deviations occur at a variety of nucleotides across all three primary cis-splicing sequences, this result is consistent with our observation that the pht1Δ strain has negative genetic interactions with genes that contribute to recognition of each of these signals.

Given the enrichment of pht1Δ-affected introns containing nonconsensus cis-splicing sequences, we wondered whether these introns were more likely to use alternative splice sites. The Pleiss group (Stepankiw et al. 2015) has recently revealed numerous instances of alternative splice site usage in S. pombe. We overlapped the list of pht1Δ-affected introns with a published list of introns using alternative 5′SSs and/or 3′SSs (Supplemental Fig. S5). Because the alternative introns were identified by lariat sequencing, the list of introns using alternative 3′SSs may also reflect alternative BPs. We found a significant enrichment of introns using alternative splice sites in the subset of pht1Δ-affected introns, the majority of which used alternative 5′SSs rather than alternative 3′SSs (BPs), suggesting that Pht1 may play an important role in 5′SS choice.

Next, we asked whether the 5′ bias of suboptimal introns affected in pht1Δ was indicative of a general distribution pattern of introns with nonconsensus cis-splicing sequences in the S. pombe genome. In order to address this question, we compared the 5′SS distance from the TSS with the 5′ SS or BP sequence strength for each intron-containing gene in the S. pombe genome. Using published position-weighted matrix scores to classify cis-splicing signal strength, we analyzed the 5′SS and BP sequence strength for all introns within 4000 base pairs (bp) of the TSS (99% of annotated introns) (Stepankiw et al. 2015). We then plotted 5′SS distances from the TSS for the subset of introns with the highest 10% and lowest 10% of 5′SS or BP scores (Fig. 3C). Importantly, introns with 5′SS or BP sequences that fell into the bottom decile for sequence strength score were more likely to be found within the first 500 bp from the TSS compared with introns with 5′SS or BP sequences in the top decile of sequence strength, implying that introns with weak cis-splicing sequences are more likely to be at the 5′ end of genes genome-wide.

Taken together, these data indicate that H2A.Z is required for efficient splicing of a subset of introns that is “weak,” characterized by suboptimal architectures, nonconsensus cis-splicing sequences, and enrichment for alternative 5′SS usage. Furthermore, introns in the 5′ ends of genes are more likely to have nonconsensus cis-splicing sequences and be retained in the absence of H2A.Z, consistent with H2A.Z's known enrichment near the TSS (Subramanian et al. 2015).

H2A.Z is enriched near introns with nonconsensus cis-splicing sequences

In order to influence splicing, one would predict that H2A.Z would be enriched near S. pombe intron–exon boundaries. To test this prediction, we used ChIP-seq (chromatin immunoprecipitation [ChIP] combined with high-throughput sequencing) to probe H2A.Z association genome-wide. We combined one replicate of our ChIP-seq with two biological replicates from a recently published data set for analysis (Clement-Ziza et al. 2014). We confirmed that the replicates were correlated highly enough to be analyzed together by plotting the average coverage of each kilobase of the genome for each data set and calculating the Pearson correlation pairwise between each replicate, showing a high correlation of all three replicates (R > 0.8) (Supplemental Fig. S6).

We first aligned all of the TSSs or intron–exon boundaries in the genome and then plotted the ChIP-seq reads that fell around each feature over the background whole-cell extract (WCE) coverage at that same position (Fig. 4A; Supplemental Fig. S7). As expected, we observed strong enrichment of H2A.Z around the TSS, with the strongest placement around the +1 nucleosome followed by the −1 nucleosome, consistent with previous observations (Zofall et al. 2009).

Figure 4.

H2A.Z is enriched at splice sites. (A) Pileup of H2A.Z ChIP-seq reads around all TSSs, 5′SSs, and 3′SSs, normalized to input sample. Graphs represent a single replicate from this study. P-values are from a K-S test for data compared with a uniform distribution. (B) Bar plots representing the average level of H2A.Z coverage measured by ChIP-seq across genes with or without introns, normalized to WCE and the coverage of both data sets. Error bars represent the SEM. (***) P < 0.001, calculated by t-test. (C) Bar plots representing the average level of H2A.Z coverage measured by ChIP-seq within ±150 bp around the TSS for intron-containing genes with or without pht1Δ-affected introns, normalized to WCE and the coverage of both data sets. Error bars represent the SEM. (***) P < 0.001 calculated by t-test. (D) Bar plots representing the average level of H2A.Z coverage measured by ChIP-seq across introns in the top and bottom deciles of 5′SS and BP strength, normalized to WCE and the coverage of both data sets. Error bars represent the SEM. (***) P < 0.001, calculated by t-test; (N.S.) not significant.

Strikingly, we also observed a strong association of H2A.Z directly upstream of both the 5′SS and 3′SS. (Fig. 4A; Supplemental Fig. S7). Similar H2A.Z association patterns for the TSS, 5′SS, and 3′SS are also observed when looking only at genes containing introns with splicing defects in the absence of pht1 at 16°C (Supplemental Fig. S8). These data suggest that H2A.Z is located in physical proximity to intron–exon boundaries in order to directly influence splicing. Moreover, the median distance from the TSS to the first intron is 249 bp in S. pombe, further correlating the location of many introns with the +1 nucleosome (Supplemental Fig. S9). However, we cannot separate the bias for H2A.Z enrichment at the +1 nucleosome and introns located near the TSS without further molecular experiments.

Because we observed an enrichment of H2A.Z near intron–exon boundaries, we next asked whether H2A.Z deposition is biased toward intron-containing genes. We calculated H2A.Z coverage measured by ChIP-seq across every gene body and compared the average levels between genes with and without introns, normalized to WCE and the average coverage of both data sets. We found significantly higher levels of H2A.Z association on intron-containing genes, providing strong correlative evidence of a functional relationship between H2A.Z and the spliceosome (Fig. 4B; Supplemental Fig. S10A).

To strengthen our hypothesis that the presence of H2A.Z promotes the splicing of proximal introns, we next analyzed whether pht1Δ-affected introns are more likely to be associated with H2A.Z. We calculated H2A.Z coverage across a 300-bp window around the TSS (150 bp upstream and downstream) and then compared the average H2A.Z association levels between genes with and without pht1Δ-affected introns, normalized as above. As expected, there is a statistically significant enrichment of H2A.Z around the TSS for introns that have splicing defects in a pht1Δ background (Fig. 4C; Supplemental Fig. S10B). Importantly, these data suggest that the presence of H2A.Z near the TSS is predictive of a weak intron and further support the claim that H2A.Z-containing nucleosomes promote splicing of nearby introns.

Given that H2A.Z is required for the efficient splicing of suboptimal introns, we next asked whether introns with weak cis-splicing sequences are more likely to coincide with H2A.Z-containing nucleosomes. We sampled introns from the top and bottom deciles of both 5′SS and BP strength scores and asked whether they had significantly different levels of normalized average H2A.Z coverage across the intron. We observed significantly higher levels of H2A.Z association on introns in the bottom decile of 5′SS strength scores compared with the top decile. However, there was not a significant difference in H2A.Z coverage between the top and bottom deciles of BP strength scores. Weak BPs lacking significant H2A.Z coverage may be due to a smaller data set of BP strength scores analyzed than of 5′SS strength scores. These data demonstrate that introns coincident with H2A.Z association are more likely to have weak 5′SSs (Fig. 4D; Supplemental Fig. S11).

Previous work in S. pombe showed a correlation between H2A.Z association and low gene expression (Zofall et al. 2009). To confirm this observation genome-wide, we compared the expression levels of H2A.Z-associated genes with the expression levels of genes without H2A.Z and found that H2A.Z-associated genes tend to be more lowly expressed (P < 0.001) (Supplemental Fig. S12). Intriguingly, it has been shown that introns in lowly expressed genes in S. pombe are more likely to undergo aberrant splicing (Stepankiw et al. 2015).

Taken together, these analyses demonstrate that H2A.Z tends to be positioned in close proximity to weak introns in lowly expressed, intron-containing genes, raising the possibility that it marks these genes to promote cotranscriptional intron recognition and splicing fidelity.

Overexpression of Prp16 suppresses H2A.Z-dependent splicing defects

Given the spatial overlap of H2A.Z and weak introns, we next wondered whether H2A.Z has any molecular connections to parts of the spliceosome involved in identifying and proofreading cis-splicing signals. Splicing ATPases known to participate in proofreading include prp28, prp11 (Prp5 in S. cerevisiae), prp16, prp22, and prp43 (Koodathingal and Staley 2013). In the S. pombe EMAP, prp28 and prp11 do not show any strong genetic interactions with the Swr1 complex, possibly due to insufficient perturbation of expression by the DAmP alleles. However, we observed synthetic sickness in double mutants of both prp22 or prp43 and pht1Δ by serial dilution growth assays at 30°C, which were exacerbated at 16°C, consistent with the strong negative interactions across the Swr1 complex with these two splicing mutants in the EMAP (Fig. 1C).

Because the synthetic sick interactions with prp43 and the Swr1 complex were particularly strong, we further characterized the prp43 allele by examining its expression level by RT-qPCR at 30°C (Supplemental Fig. S13A). prp43 is overexpressed in the prp43 strain by more than three times the level of wild type. As Prp43 promotes the discard of suboptimal substrates, the synthetic sickness exacerbated at 16°C in the double-mutant pht1Δ prp43 strain might be explained by increased spurious discard of weak introns (Mayas et al. 2010).

Because the prp16 strain was not present in the S. pombe EMAP analysis, we generated double mutants with the pht1Δ or swr1Δ strain by hand and performed serial dilution growth assays across different temperatures to look for genetic interactions (Fig. 5A). Prp16 is a splicing DEAH-box helicase required for transition from the first catalytic step to the second catalytic step of splicing and shown in vitro to possess proofreading activity, acting at the 5′SS and BP (Horowitz 2011). Strikingly, the pht1Δ prp16 and swr1Δ prp16 double-mutant strains suppressed both the pht1Δ cold sensitivity and the swr1Δ heat sensitivity. RT-qPCR revealed that prp16 transcript levels are at wild-type level in pht1Δ but more than threefold above the wild-type level in prp16 and pht1Δ prp16 strains at 30°C and 16°C (Supplemental Fig. S13A,B).

Figure 5.

Overexpression of prp16 ATPase suppresses growth and splicing defects in pht1Δ and swr1Δ at restrictive temperatures. (A) Serial dilutions (1:5) of wild-type and single- and double-mutant strains grown at 16°C, 30°C, and 37°C. (B) Heat maps of the splicing profiles of the pht1Δ or swr1Δ, prp16, and their respective double-mutant strains. Cultures were grown to mid-log phase at 30°C and shifted for 9 h to 16°C or for 2 h to 37°C where indicated; cDNA from each single- and double-mutant strain was competitively hybridized on the splicing-specific microarray against that from wild type. The stacked bar plot shows the number of introns with a splicing index score >0.3 or less than −0.3. The color scale shows distribution of the severity of the splicing defect. (N.S.) Not significant

Because overexpression of prp16 suppresses both pht1Δ and swr1Δ temperature-dependent growth defects, we next tested whether the double mutants also suppress the intron accumulation defects of pht1Δ and swr1Δ at 16°C and 37°C, respectively. We used the splicing-specific microarray platform to assess the splicing efficiency of prp16 and pht1Δ prp16 at 16°C and prp16 and swr1Δ prp16 at 37°C (Fig. 5B). The prp16 strain by itself had a mild splicing defect at 16°C with 190 introns affected and at 37°C with 83 introns affected, which is not as strong as either pht1Δ or swr1Δ, respectively. However, overexpression of prp16 in strains lacking pht1 or swr1 led to a significant reduction in the total number of introns with splicing defects. In pht1Δ prp16, there was a partial reduction in the total number of introns affected, from 883 in pht1Δ to 500, consistent with the partial suppression of the cold-sensitive growth defect.

Importantly, the vast majority of the subset of 883 introns retained in pht1Δ at 16°C is rescued in the prp16 strain (97%) and the pht1Δ prp16 strain (85%) at 16°C. We further analyzed the subset of introns with significant splicing changes in strains overexpressing prp16 at 16°C and found 167 introns in the prp16 strain and 370 introns in the pht1Δ prp16 strain that are distinct from the introns that accumulate in the absence of pht1. This subset of prp16-affected introns is too small to perform intron feature analysis on with the method that we used. These data suggest that excess Prp16 is not favorable for constitutive splicing of all introns. In the swr1Δ prp16 strain, we detected no splicing events with significant changes from wild type, indicating that overexpression of prp16 restored splicing of all 155 introns that accumulated in swr1Δ at 37°C. However, because of improved growth in the double mutants, we cannot rule out the possibility of indirect suppression of intron accumulation.

Because the splicing of weak introns in pht1Δ at 16°C can be rescued by overexpression of the proofreading ATPase Prp16, we propose that nonconsensus introns require additional inspection by the spliceosome to be properly processed in suboptimal environmental conditions. Furthermore, because overexpression of Prp43 exacerbates the pht1Δ-dependent splicing defect, these data suggest that H2A.Z may have a role in preventing spurious discard of suboptimal substrates by promoting the splicing reaction.

Restoration of consensus cis-splicing sequences suppresses the H2A.Z-dependent splicing defect

Next, we sought to directly test the requirement for H2A.Z in promoting splicing of “weak” introns. If H2A.Z indeed promotes splicing of introns with nonconsensus cis-splicing sequences, then we would predict it to be dispensable for splicing when the intron's cis-splicing sequences are strengthened to consensus.

To address this question, we chose a model intron from the single-intron-containing gene SPBPJ4664.05, which is found at the 5′ end of the gene and accumulates in the absence of H2A.Z by microarray (splicing index = 1.187; P < 0.0006). SPBPJ4664.05 is enriched for H2A.Z deposition by ChIP-seq and is lowly expressed (Marguerat et al. 2012). SPBPJ4664.05 was tagged with an auxotrophic marker, and the nonconsensus 5′SS and BP sequences of its intron were mutated to the S. pombe consensus sequences (Fig. 6A). These constructs were integrated into the genome, replacing the endogenous alleles in wild-type and pht1Δ strain backgrounds. Growth conditions were consistent with the 16°C microarray conditions, and RT-qPCR was performed to assess splicing efficiency of SPBPJ4664.05 intron 1.

Figure 6.

Mutating cis-splicing sequences to consensus suppresses the splicing defect of a model intron in pht1Δ at 16°C. (A) Cartoon depicting the architecture and 5′SS and BP sequences of intron 1 of model gene SPBPJ4664.05. Mutation of cis-splicing sequence nucleotides to the consensus sequence is indicated in red. (B) Bar graphs representing the splicing efficiency of SPBPJ4664.05 intron 1 in wild-type and pht1Δ strains with and without mutations to the 5′SSs and BPs as indicated above. Cultures were grown to mid-log phase at 30°C and shifted for 9 h to 16°C. Bars represent the ratio of intronic to exonic signal measured by RT-qPCR. Error bars represent the SEM of three biological replicates. (***) P < 0.001, calculated by one-tailed t-test.

Importantly, we found that strengthening the cis-splicing sequences of SPBPJ4664.05 intron 1 to consensus in the pht1Δ background is sufficient to restore splicing efficiency to the wild-type level (Fig. 6B). Together with our genome-wide microarray data and global H2A.Z ChIP, these results strongly suggest that H2A.Z promotes the efficient removal of weak introns by the spliceosome.

Discussion

Proper recognition and removal of introns requires regulation at a number of steps. How cells integrate multiple inputs to make cotranscriptional splicing decisions is still an area of active investigation. In this study, we demonstrated a requirement for the histone variant H2A.Z and the Swr1 nucleosome remodeling complex in the efficient splicing of weak introns. We showed that H2A.Z is enriched at intron–exon boundaries near promoters of lowly expressed genes. Additionally, introns located at the 5′ end of S. pombe genes are more likely to contain nonconsensus cis-splicing signals. Correlation between H2A.Z occupancy and nonconsensus introns suggests that the histone variant itself may mark introns with weak splicing signals to aid in their recognition and removal. Consistent with this idea, we show that overexpression of a proofreading ATPase, Prp16, rescues growth and splicing defects in the absence of H2A.Z-containing nucleosomes. We propose a model in which the presence of H2A.Z proximal to introns alters the kinetics of cotranscriptional intron recognition and proofreading, thereby promoting the splicing of introns with weak splice sites.

H2A.Z promotes the recognition and removal of weak promoter-proximal introns

The genome of S. pombe is rich with thousands of cis-splicing signals that are not constrained to a consensus sequence (Wood et al. 2002; Kuhn and Kaufer 2003). Such diversity in sequences implies a need for multiple layers of regulation for every splicing decision. The work presented here strongly implicates H2A.Z and the Swr1 complex in facilitating these decisions in fission yeast. Our results are consistent with published studies from the Pleiss group (Larson et al. 2016) in S. pombe and the Stevens group (Sorenson and Stevens 2014) in S. cerevisiae, in which multiple components of the Swr1 complex were identified in genome-wide screens for splicing efficiency changes of a model gene construct.

One way that the Swr1 complex may promote splicing is through the recruitment of splicing factors to the nascent pre-mRNA during spliceosome assembly and/or catalysis. Our work shows that H2A.Z, in addition to being enriched around the TSS, is also more likely to be found in nucleosomes flanking weak introns in lowly expressed genes, spatially correlating H2A.Z deposition and splicing at genes with weak gene expression signals. In humans, SF3B was shown to weakly coimmunoprecipitate with H2A.Z, suggesting that there may be physical interactions underlying these genetic observations (Tolstorukov et al. 2012). In addition, mass spectrometry experiments of the Swr1 complex in S. pombe and S. cerevisiae show that other early splicing factors are present at subcomplex levels, particularly splicing factors for which there are strong negative genetic interactions in both the S. pombe and S. cerevisiae EMAPs (Shevchenko et al. 2008; Wilmes et al. 2008; Kim et al. 2009). Further biochemistry is needed to examine the nature of these genetic interactions.

Alternatively, and not mutually exclusively, H2A.Z may help coordinate splicing and transcription by influencing the kinetics of RNAP II initiation and elongation through promoter-proximal introns. To further understand H2A.Z's role in the cotranscriptional recruitment of the spliceosome and/or dynamics of transcription, we performed ChIP-qPCR with splicing factors (U1 and U2AF) and RNAP II (holoenzyme and Ser5-P) at multiple genes of interest that we observed to have both H2A.Z deposition and intron retention in pht1Δ at 16°C (data not shown). However, because of the correlation between H2A.Z occupancy and low levels of gene expression, we could not consistently detect either spliceosome proteins or RNAP II at levels over background at genes with pht1Δ-affected introns. These results are consistent with previous reports in S. cerevisiae that show that RNAP II has a low ChIP signal in the htz1Δ (H2A.Z) background (Adam et al. 2001).

Recent work from the Johnson group (Neves et al. 2017) in S. cerevisiae shows that, in the absence of H2A.Z, the Ser2-P population of RNAP II, representing elongating RNAP II, accumulates at the 5′ ends of the gene bodies. They also observed significant perturbations in snRNP association at introns with nonconsensus BPs, suggesting a defect in spliceosome assembly, rearrangements, or discard. Taken together with our data, these observations demonstrate that H2A.Z's influence on transcription is correlated with its role in splicing of nonconsensus introns, possibly through changes in the kinetics of spliceosome rearrangements. Observations made regarding RNAP II elongation in S. cerevisiae are likely applicable to the fission yeast system. Previous analyses of our EMAPs showed that Swr1 complex components as a group had highly correlated global genetic interaction profiles in both S. pombe and S. cereviaise, strongly suggesting that this protein complex plays similar roles in the two yeast species.

H2A.Z likely influences splicing fidelity

The splicing reaction proceeds in a complex stepwise manner, driven forward by ATPases that regulate the rate of spliceosomal rearrangements, allowing proofreading of the substrates (Staley and Guthrie 1998; Koodathingal and Staley 2013). Given our observation that weak, nonconsensus introns are sensitive to loss of H2A.Z, we were curious as to whether expression of proofreading ATPases—the bastions of splicing fidelity—could compensate for loss of H2A.Z at introns with weak and/or alternative cis-splicing signals. Indeed, we observed strong suppression of the pht1Δ-dependent splicing defect in cells overexpressing the splicing fidelity factor DEAH-box ATPase Prp16. This requirement for additional Prp16 is underscored in our experiments performed at 16°C, where the loss of pht1 likely contributes to slow/inefficient RNAP II elongation that is further exacerbated at low temperature. Under these conditions, overexpression of Prp16 may rescue splicing by increasing the association frequency of Prp16 with the spliceosome, thereby increasing the probability that Prp16 can sample and engage the cis-splicing signals of weak, promoter-proximal introns (Fig. 7). By increasing the likelihood that Prp16 is available through an increase in its local concentration, both first- and second-step catalysis can be pushed forward even in the absence of H2A.Z to mark weak introns and coordinate cotranscriptional intron recognition and removal. This role for Prp16 is consistent with a report from the Cheng group (Tseng et al. 2011), which showed that Prp16 has an ATP-independent role in promoting the first catalytic step of splicing, and multiple reports from the Staley group (Koodathingal et al. 2010), which characterized Prp16's roles in kinetically proofreading the second step, and allows for sampling of nonconsensus splice sites by the spliceosome (Semlow et al. 2016)

Figure 7.

Model illustrating the effect of H2A.Z loss on splicing efficiency and suppression by overexpression of prp16 and exacerbation by overexpression of prp43. Cartoon depicting a gene with a weak intron near the +1 nucleosome. When H2A.Z is absent, splicing of this weak intron does not proceed efficiently. Overexpression of Prp16 ATPase recues this splicing efficiency defect, possibly through altering the kinetics of splicing catalysis. Through mass action, overexpressed Prp16 is more likely to be bound to the spliceosome, thereby preventing discard by Prp43, promoting weak splice site sampling, and/or driving splicing catalysis forward. Overexpression of Prp43, on the other hand, preferentially acts to discard the spliceosomes on the weak substrate, which are possibly stalled due to the absence of H2A.Z.

Such a model predicts that, in a wild-type cell, splicing of suboptimal substrates requires additional attention from the spliceosome, which is afforded through coordination between H2A.Z and the spliceosome. By allowing the spliceosome to try and fail several times on sub-optimal substrates before the first catalytic step is completed, the cell can avoid discarding substrates with weak splice sites. Consistent with this idea, overexpression of Prp43 in the pht1Δ strain causes synthetic sickness, possibly through increased association of Prp43, leading to discard of already poorly spliced, weak introns. Rescuing the growth and splicing defects of pht1Δ by overexpressing Prp16 suggests that H2A.Z has a role in promoting the splicing of weak introns by prompting additional inspection by the spliceosome, thus preventing spurious discard of nonconsensus splice site-containing introns. This model is bolstered by our observation that mutating an intron with nonconsensus 5′SS and BP sequences to consensus suppresses the requirement for H2A.Z in promoting efficient splicing.

Taken together with results from the Johnson group (Neves et al. 2017), our genetic and molecular data implicate the H2A.Z histone variant as a signal that marks weak introns to promote cotranscriptional splicing. By uncovering the complexities of how chromatin influences splicing fidelity and efficiency, we can begin to interrogate how uncoupling these coordinated processes may result in cellular dysfunction and disease.

Materials and methods

Strains

A list of strains used in this study is available in Supplemental Table S3.

EMAP analysis

EMAP data and their generation were as described previously (Patrick et al. 2015; Ryan et al. 2012). Evaluation of the links between the Swr1 complex and GO processes or complexes was performed as in Patrick et al. (2015) using the same data set. Additional complexes related to splicing were defined by manual curation and are available in Supplemental Table S1.

RNA isolation and qRT–PCR

Strains were grown in YES5 overnight at 30°C and then diluted to 0.1 OD600 the next morning to allow for growth to mid-log phase. Temperature shifts occurred when strains had reached 0.3–0.5 OD600. Shifts to 37°C took place for 2 h, and shifts to 16°C took place for 9 h. Cultures were harvested by centrifugation between 0.5 and 0.8 OD600. Pellets were washed once in water and then flash-frozen in liquid nitrogen and stored at −80°C. RNA was isolated using hot acid phenol followed by isopropanol precipitation as described previously (Bergkessel et al. 2011).

For RT-qPCR, RNA was treated with DNase I (Fermentas) and then MMLV1-RT (Promega) with 0.5 mg/mL dN9 primers to generate cDNA. The remaining RNA was then removed from cDNA by NaOH treatment, and the cDNA was purified using a spin kit (Zymo). qPCR was performed using a CFX96 real-time system (Bio-Rad).

Splicing-specific microarray data collection and feature analysis

Collection and analysis of splicing-specific microarray data was performed as described previously (Lipp et al. 2015; Patrick et al. 2015). Logistic regression analysis and subsequent intron feature analysis were performed as described previously (Lipp et al. 2015). Microarray data for this study can be accessed at the NCBI Gene Expression Omnibus under accession number GSE97984.

ChIP-seq library construction

ChIP was performed as described previously (Patrick et al. 2015). After a ChIP reaction, the sequencing library preparation was performed as described previously (Dumesic et al. 2015) with the following modifications: The DNA yield of a single ChIP reaction—or 100 ng of WCE DNA—was used as starting material. Library concentration was determined by Qubit (Invitrogen).

Sequencing data processing, read alignment, and analysis

Alignment of reads, generation of BAM files, and generation of Bedgraph files were performed as described previously (Dumesic et al. 2015). Bedgraph files were used for analysis of signal around intron–exon junctions as described below. Bedgraph files were visualized using the Integrative Genomics Viewer 2.0.30 (Broad Institute). The S. pombe genome from Pombase (http://www.pombase.org) was used as the reference.

Additional ChIP-seq data of wild-type and pht1-3xFlag samples were obtained from ArrayExpress with identifier E-MTAB-2650 (Clement-Ziza et al. 2014). These data were treated as additional replicates and processed in the same manner as the library generated in this study. To confirm that replicates were similar enough to be used together in analysis, the average coverage of each kilobase of the genome for each data set was plotted. The Pearson correlation was calculated for each replicate compared with another replicate using these average kilobase coverage values.

ChIP-seq data for this study can be accessed at the NCBI Gene Expression Omnibus under accession number GSE97984.

Meta-intron–exon boundary and meta-TSS graphs

Each intron–exon boundary—either 5′SS, 3′SS, or both—for every intron in the genome as well as the subset of introns with splicing defects in the pht1Δ microarray at 16°C was aligned, and the pileup of ChIP signal 500 bp upstream of and downstream from this feature was analyzed. The total coverage around this type of feature was then normalized on a per-base basis to that of the total coverage of the input sample at these same locations. To control for differences in sequencing depth, each sample was further normalized to read depth such that it equaled read depth in the input sample. Meta-TSS plots were generated in a similar manner by aligning the TSS for each gene in the genome and analyzing ChIP signal 500 bp upstream of and downstream from the feature, normalized to input and coverage of both data sets. The significance of ChIP enrichment in these pileup graphs was determined using the Kolmogorov-Smirnov (K-S) test compared with input sample signal distribution.

ChIP-seq coverage analysis

For each gene in the genome, the number of reads in H2A.Z ChIP-seq data sets aligning to the gene body or the intron were counted. The same process was applied to the input sample data set. Fold enrichment was calculated for each gene by dividing ChIP counts by input data set counts, normalized to the total read counts for both ChIP and input samples: (ChIP data set reads for gene/input data set reads for gene) × (input data set total reads/ChIP data set total reads). The same analysis was applied for each intron and TSS region in the genome, with TSS region defined as 150 bp upstream of and downstream from each TSS. Mean enrichment of ChIP signal over input was compared between genes containing introns whose splicing was altered by the presence/absence of H2A.Z, between introns that are sensitive to the presence/absence of H2A.Z, or between the TSS region of genes that contain an intron sensitive to the presence/absence of H2A.Z. This analysis was repeated on introns that were found to be in the top and bottom deciles of 5′SS and BP scores in the S. pombe, as reported previously (Stepankiw et al. 2015).

Analysis of intron position and strength in the S. pombe genome

Intron locations were obtained from Pombase (http://www.pombase.org). S. pombe 5′SS and BP strength scores were reported previously (Stepankiw et al. 2015). These data consisted of position-weight matrix scores for the majority of annotated 5′SSs (5294) and BPs (3076). Some introns do not have scores for these cis-splicing sequences because the methods used were not able to positively identify these features. The distance of each intron from the TSS of its gene was graphed on a theoretical gene body. Distributions across a theoretical gene body were also graphed for introns in the top and bottom deciles of 5′SS and BP strength scores. The significance of intron distribution difference was determined by the K-S test.