Position-specific intron retention is mediated by the histone methyltransferase SDG725
BMC Biology volume 16, Article number: 44 (2018)
Intron retention (IR), the most prevalent alternative splicing form in plants, plays a critical role in gene expression during plant development and stress response. However, the molecular mechanisms underlying IR regulation remain largely unknown.
Knockdown of SDG725, a histone H3 lysine 36 (H3K36)-specific methyltransferase in rice, leads to alterations of IR in more than 4700 genes. Surprisingly, IR events are globally increased at the 5′ region but decreased at the 3′ region of the gene body in the SDG725-knockdown mutant. Chromatin immunoprecipitation sequencing analyses reveal that SDG725 depletion results in a genome-wide increase of the H3K36 mono-methylation (H3K36me1) but, unexpectedly, promoter-proximal shifts of H3K36 di- and tri-methylation (H3K36me2 and H3K36me3). Consistent with the results in animals, the levels of H3K36me1/me2/me3 in rice positively correlate with gene expression levels, whereas shifts of H3K36me2/me3 coincide with position-specific alterations of IR. We find that either H3K36me2 or H3K36me3 alone contributes to the positional change of IR caused by SDG725 knockdown, although IR shift is more significant when both H3K36me2 and H3K36me3 modifications are simultaneously shifted.
Our results revealed that SDG725 modulates IR in a position-specific manner, indicating that H3K36 methylation plays a role in RNA splicing, probably by marking the retained introns in plants.
Intron retention (IR), a specific form of pre-messenger RNA (pre-mRNA) alternative splicing (AS), has attracted increasing attention given its role in global gene expression regulation in both animals and plants [1,2,3,4,5,6]. Notably, IR is the most prevalent form of AS in higher plants, possibly due to the shorter intron length in plants than in animals [7, 8]. Genome-wide analyses revealed that greater than half of the AS events in rice belong to IR [9, 10].
Several studies have highlighted the functional importance of IR in plants. For example, IR has been shown to associate with abiotic stress response in barley [5, 11]. In strawberry, the level of IR is significantly reduced in post-fertilization compared to pre-fertilization, suggesting the involvement of IR in fruit maturation . Retention of an intron in the 5′ UTR of the Zinc-Induced Facilitator 2 gene (ZIF2) enhances zinc tolerance in Arabidopsis . Two intron-retained transcripts in Arabidopsis have been shown to remain in the nucleus to avoid nonsense-mediated degradation (NMD) . During Arabidopsis gametophyte development, IR regulates translation in a transcription-independent and spliceosome-dependent manner . Taken together, all these findings underscore the importance of IR in plant growth and development.
Epigenetic regulators are known to be involved in AS regulation [4, 15, 16]. Histone modifications are particularly interesting because of their potential links between chromatin structure and co-transcriptional pre-mRNA splicing. Chromatin immunoprecipitation sequencing (ChIP-seq) analyses in human and mouse cells showed that H3K36me3 (tri-methylation of histone H3 at lysine 36) signals in introns and alternatively spliced exons are considerably lower than those in constitutive exons, suggesting that H3K36me3 modification likely acts as a mark for exons in animal cells . Furthermore, H3K36me3 participates in AS by recruiting the splicing factor polypyrimidine tract-binding protein (PTB) via MRG15, an H3K36me3 reader protein in human . BS69, a specific reader protein for H3.3K36me3, is involved in pre-mRNA splicing, especially IR, by interaction with the U5 small cytoplasmic fractionation extraction ribonucleoprotein (snRNP) in human cells . Pajoro et al. discovered that H3K36me3 played a role in AS and flowering control in Arabidopsis, wherein mutants of SDG8 and SGD26, two methyltransferases of H3K36, affect temperature-dependent flowering . All these findings support the notion that H3K36me3 plays a direct role in regulating AS.
However, it remains unclear whether all three forms of H3K36 methylation (mono-, di-, and tri-) are involved in AS regulation, especially IR, in plants. We previously reported that two H3K36-specific methyltransferases, SDG725 and SDG708, modulate gene transcription and affect rice growth and development [19,20,21]. In this study, we investigated splicing alterations in the SDG725- and SDG708-knockdown rice mutants by RNA sequencing (RNA-seq). We found that knockdown of SDG725 led to altered IR in thousands of genes. In addition, IR events tend to increase at the 5′ portion but decrease at the 3′ part of the gene body when comparing the SDG725-knockdown mutant to the wild-type (WT) plants. The results coincided with a higher H3K36me2 occupancy at the 5′ part but a lower one at the 3′ part of the gene body. In contrast, SDG708 knockdown did not cause either these promoter-proximal shifts of IR or histone modification. Our work discovered a previously unknown shift of IR and its possible epigenetic regulator in rice.
SDG725 regulates a global shift of intron retention in rice
Since H3K36 methylation has been proposed for splicing regulation in animals, we sought to investigate whether a similar mechanism is also employed in plants. We took advantage of two transgenic rice lines previously generated, in which SDG725 or SDG708 was efficiently knocked down by RNA interference [19, 21]. Quantitative mass spectrometry showed that knockdown of SDG725 led to an increased level of H3K36me1 modification, but decreased levels of both H3K36me2 and H3K36me3 (Additional file 1: Figure S1) . Two biological replicates of RNA-seq libraries were constructed, sequenced, and analyzed for 725Ri-1 (a stable RNAi line of SDG725), 708Ri-1 (a stable RNAi line of SDG708), and WT rice plants (Additional file 2: Table S1) [19, 21]. As the result of SDG725 knockdown, RNA-seq analyses revealed that 462 and 496 genes were up- and down-regulated, respectively (Additional file 1: Figure S2a). Gene ontology analysis showed that the differentially expressed genes (DEGs) in 725Ri-1 are enriched in metabolic and biosynthetic processes (Additional file 1: Figure S2b). The DEGs in 708Ri-1 (245 up- and 222 down-regulated) also are enriched in metabolic processes, but only a small fraction of them overlapped with those found in 725Ri-1 mutant plants , indicating the distinct biological roles of these two H3K36-specific methyltransferases in rice.
The splicing changes in rice were further analyzed by the SplAdder approach . Consistent with previous findings , IR is the predominant form of AS events identified by comparing to rice genome annotation (Additional file 2: Table S2). To perform a quantitative investigation on IR, we used the intron retention index (IRI) to estimate the extent of retention for each annotated intron . We found that 4714 genes contain one or more altered IR, including 2089 IRI-up introns (≥ twofold increase in IRI) and 4214 IRI-down introns (≥ twofold decrease in IRI) by comparing the 725Ri-1 plants to the WT rice plants. Ten IRI-altered introns were selected for validation by reverse transcription followed by quantitative polymerase chain reaction (qRT-PCR). Primer pairs were designed to detect either spliced or intron-retained transcripts for each IR event (Additional file 2: Table S3). For all 10 cases, qRT-PCR results confirmed the RNA-seq data (Fig. 1a). To our surprise, we found that the IRI-up introns were favored at the 5′ portion of the gene body while the IRI-down introns were preferred at the 3′ portion of the gene body when comparing 725Ri-1 to WT rice plants using stringent coverage requirements (Fig. 1b, Additional file 1: Figure S3). The conclusion remained the same when the requirements of intron length coverage were relaxed (Additional file 1: Figure S4). In contrast, no location bias was observed between the 3326 IRI-up and 2497 IRI-down introns in 708Ri-1 (Fig. 1c, Additional file 2: Table S4). For exon skipping, no obvious position-specific splicing alteration was found by comparing 725Ri-1 with WT rice. Our results showed that knockdown of SDG725 but not SDG708 alters IR in a position-specific manner.
Genes with increased IR show reduced expression levels
Among all the genes with one or more retained introns, 1399 genes only contain IRI-up introns, 2908 genes only contain IRI-down introns, and 407 genes contain both IRI-up and IRI-down introns (Fig. 2a). IR may regulate gene expression through several different mechanisms, including but not limited to the following: (1) intron-retained transcripts are unstable, being degraded by nuclear RNA surveillance machinery or cytoplasmic NMD [3, 13, 24,25,26]; (2) intron-retained transcripts are stable and can be translated into new protein variants, where those introns are known as “exitrons” [27,28,29]; (3) a retained intron in a 5′ UTR regulates translation initiation  while IR in a 3′ UTR can either repress mRNA stability and translation by introducing more regulatory cis elements  or stabilize mRNA by avoiding NMD . Therefore, global IR coupled with RNA stability could serve as a novel mechanism to fine-tune gene expression in rice. RNA-seq data were then used to determine the effect of altered IR on gene expression. We found that the expression level of the genes with only IRI-down introns tends to be up-regulated compared to those with only IRI-up introns, and the expression changes of the 407 genes with both IRI-up and IRI-down introns fall in between (Fig. 2b). These observations agreed with the notion that IR may serve as a post-transcriptional mechanism to reduce gene expression [1,2,3,4]. Consistent with Fig. 1b, in the 407 genes, IRI-up introns show a similar preference at the 5′ part of the gene body and IRI-down introns are enriched at the 3′ end of the gene body (Fig. 2c). The same occurs with the genes with IRI-up only and IRI-down only introns (Fig. 2d).
Intron retention shift correlates with distribution shifts of H3K36 methylations
We next asked why IR is shifted as the result of SDG725 knockdown. Considering the role of H3K36 methylation in IR in animals, we acquired H3K36 mono-, di-, and tri-methylation (H3K36me1/me2/me3) profiles of the 725Ri-1, 708Ri-1, and WT rice plants using ChIP-seq (Additional file 2: Table S5). By looking at gene-dense regions we found that H3K36me2 and H3K36me3 profiles exhibit apparent promoter-proximal shifts in 725Ri-1 but not in 708Ri-1 compared to those in WT rice plants (Fig. 3a). These observations were further confirmed by genome-wide analyses (Fig. 3b–d). It is worth noting that H3K36me1 levels are increased in nearly the entire gene body (Fig. 3b, left panel) upon SDG725 knockdown. In contrast, SDG708 knockdown led to a global decrease of all three H3K36 methylation marks (Fig. 3a, b).
To obtain a detailed view of H3K36 methylation changes at the individual gene level, we computed the difference of ChIP-seq tag density gene by gene between knockdown and WT plants (see details in Methods, Additional file 1: Figure S5). For almost all transcribed loci in the 725Ri-1 plants, H3K36me1 levels are increased across the gene body (Fig. 3c, left panel), while H3K36me2 (Fig. 3c, middle panel) and H3K36me3 (Fig. 3c, right panel) levels are increased at the 5′ end but decreased at the 3′ end of the gene body. However, this histone methylation positional bias was not detected in the 708Ri-1 plants (Fig. 3d), suggesting that the promoter-proximal shifts of H3K36me2/me3 are specific for the 725Ri-1 plants.
H3K36me2/me3 shifts positively correlate with IR shifts caused by SDG725 knockdown
SDG725 knockdown caused global shifts of H3K36me2 and H3K36me3, providing us with a unique opportunity to examine the potential contributor for IR shift. We therefore divided the transcribed genes into four groups. For type I genes, only H3K36me2 show a promoter-proximal shift. Type II genes are designated for H3K36me3 shift towards the 5′ end, while type III genes exhibit pattern shifts for both H3K36 methyl marks. Type IV genes show no shift for either H3K36me2 or H3K36me3. We noticed that H3K36me2 shift alone led to a considerable effect on the location bias of IR events (Fig. 4a). Similarly, H3K36me3 shift alone has an impact on IR shift (Fig. 4b). Intriguingly, when both H3K36me2 and H3K36me3 profiles are simultaneously shifted, IR showed a more significant shift compared to either H3K36me2 or H3K36me3 alone (Fig. 4c). As a control, unshifted H3K36me2/me3 exhibit no significant effect on IR switch (Fig. 4d). These results suggested that both H3K36me2 and H3K36me3 are involved in modulating IR shift, although a collaborative mechanism may exist in regulating IR.
However, the signals of H3K36 mono-, di-, and tri-methylation were also compared in differentially retained introns. Interestingly, both the levels of H3K36me2 and H3K36me3 at IRI-up introns were higher than those at IRI-down introns (Additional file 1: Figure S6), supporting the notion that H3K36me2/me3 probably demarcate IR in rice.
Validation of the association between intron retention and H3K36me2 modification
The candidate gene approach was used to validate the potential link between H3K36me2 and IR. Two genes (LOC_Os01g24680, encodes putative 3-hydroxyacyl-CoA dehydrogenase; LOC_Os08g01760, encodes putative nutrition dehydrogenase) were selected from the 407 genes containing both IRI-up and IRI-down introns. For each gene, we selected six regions along the gene body to test H3K36me2 occupancy by ChIP-PCR. In addition, two IR events (one IRI-up and one IRI-down) for both genes were also examined by qRT-PCR in the same samples by multiple primer pairs spanning both intronic location and donor/acceptor sites (Fig. 5a, b, Additional file 2: Table S6). The results confirmed that the chromatin regions with increased levels of H3K36me2 modification produce increased IR events; accordingly, the chromatin regions with reduced H3K36me2 occupancy are associated with decreased IR events (Fig. 5).
Transcripts with changed IR mainly accumulate in the nucleus
Cellular fractionation was performed to examine whether the IR transcripts are accumulated in the nucleus or cytoplasm. The success of cellular fractionation was confirmed by qRT-PCR of genes specifically expressed in the nucleus or cytoplasm (Additional file 1: Figure S7) [33, 34]. RNA-seq analysis on each fraction showed that transcripts with changed IR in rice dominantly accumulated in the nucleus (Fig. 6a, b, Additional file 2: Table S7). To further examine whether cytoplasmic NMD contributes to the transcript abundance, both 725Ri-1 and WT plants are treated with cycloheximide (CHX) to stabilize transcripts that are otherwise degraded by NMD, followed by RNA-seq analysis. Transcripts with a premature termination codon (PTC) showed an increased steady-state expression level upon CHX treatment (Additional file 1: Figure S8), indicating the successful inhibition of NMD. In contrast, genes with IR did not show an obvious expression increase upon CHX treatment (Additional file 1: Figure S9), suggesting cytoplasmic NMD has limited contribution to overall expression of IR loci in rice. These results indicate that intron-retained transcripts in rice are mainly accumulated in the nucleus rather than the cytoplasm and contribute to steady-state gene expression.
IR is thought to play a critical role in gene expression regulation in animals and plants. However, how IR is regulated in plants remains largely unknown. In this study we revealed that SDG725 may regulate global IR through H3K36me2/me3 modifications. We previously reported that SDG725 acts as an H3K36 methyltransferase and functions in promoting gene transcription [19, 20]. In both WT and 725Ri-1 plants, a positive correlation was observed between transcript abundance and H3K36me1/me2/me3 levels (Additional file 1: Figure S10a), extending previous findings obtained from other species [17, 35]. We discovered that changes in average H3K36me2/me3 occupancy level positively correlate with expression changes between 725Ri-1 with WT plants (Additional file 1: Figure S10b, c, left panels). The shift of H3K36me2/me3 affects IR but not expression level if the overall H3K36me2/me3 level does not change (Additional file 1: Figure S10b, c, right panels; Additional file 1: Figure S11). Taken together, we propose that SDG725 may control gene expression through two different mechanisms: (1) by modulating gene transcription via changing the overall levels of H3K36me2/me3, and (2) by regulating IR shift through H3K36me2/me3 shifts. How H3K36me2 or H3K36me2/me3 modulates IR in rice remains an open question. It might be achieved by reducing Pol II elongation rate and/or through a chromatin adaptor mechanism. In the first model, increased H3K36me2/me3 levels may slow down the elongation of Pol II. A longer dwell time of Pol II on introns may recruit splicing repressive factors or inhibit positive splicing factors to promote IR . Alternatively, H3K36me2/me3 may be recognized by a specific “reader” protein, which interacts with splicing repressive factors to promote IR [4, 16]. Notably, the two models are not mutually exclusive and may act in concert to recruit splicing regulators. While significant expression change was not detected for known splicing factors between 725Ri-1 and WT rice (Additional file 2: Table S8), further investigations are warranted to identify H3K36me2/H3K36me3 reader(s) as well as downstream factors functionally involved in splicing regulation in rice.
The coupling of transcription and splicing is prevalent . We observed a positive correlation between efficient splicing (or less IR) and gene expression level (Fig. 2b), indicating that IR negatively impacts the steady-state RNA level likely due to a higher degradation of IR transcripts. Moreover, increased IR events in the 5′ half of the gene in the SDG725-knockdown line suggest interaction between transcription and splicing machineries. Besides IR, other factors such as transcription activity are also possible contributors to steady-state mRNA expression.
To investigate the characteristics of retained introns, we examined several features including intron length, GC content, and splice site strength. Compared with spliced introns, retained introns tend to have much longer length, higher GC content, and weaker splice strength in the present rice study (Additional file 1: Figures S12–S14) [2, 3, 37]. Notably, the retained introns in animals tend to be shorter compared to spliced introns . Because the average intron size is much bigger in animals (9519 bp in mouse and 11,538 bp in human) than in plants (407 bp for rice), we speculate that in animals a shorter intron is more likely to be recognized as an exon and has a higher tendency to be retained, whereas in plants a longer intron has a higher tendency to be recognized as an exon and is subsequently retained. The underlying mechanism is expected to be complicated and deserves further characterization.
The question of why the SDG708 knockdown shows a global decrease of H3K36me2/me3 while the SDG725 mutant displays a pattern shift is intriguing. One possible explanation is the different enzyme specificity between SDG725 and SDG708. Although both are H3K36 methyltransferases, SDG725 plays a major role in mono- to di- and di- to tri-methylation, while SDG708 functions more on me0 to mono-methylation according to our previous studies [19, 38]. As expected, knockdown of SDG708 reduced the level of H3K36me1 around the transcription start site (TSS) and transcription termination site (TTS) (Fig. 3b). The occupancy of H3K36me2 and H3K36me3 was also reduced (Fig. 3b), possibly due to the lack of H3K36me1. In addition, SDG725 but not SDG708 contains a CW domain, which can bind to H3K4 methylation that is enriched at promoter regions . When SDG725 became limiting under the knockdown condition, it was preferentially recruited to the 5′ end of the gene, thereby showing more reductions of di- and tri-methylation of H3K36 at the 3′ portion of the gene body. This explains that the ChIP-seq result, represented as a relative distribution instead of an absolute level, showed a relative decrease of H3K36me2/me3 level in the 3′ half of the gene body but a relative increase of H3K36me2/me3 level close to the 5′ end (Fig. 3b).
Several lines of evidence support the notion that H3K36me2 in rice probably serves as the functional counterpart of H3K36me3 in animal cells [4, 16, 17]. H3K36me2 in plants and H3K36me3 in animals share similar occupancy profiles at transcribed regions (Fig. 3) . The H3K36me3 patterns in both rice and Arabidopsis resemble the H3K4me3 or H3K79me3 profile in animals (Additional file 1: Figure S15) [40,41,42,43]. These observations suggest that plants and animals may employ distinct epigenetic factors (e.g., writers and reads) to coordinate transcription and post-transcriptional gene regulation, although the underlying regulatory principles are conserved during evolution. Our study demonstrated for the first time that a histone methyltransferase can regulate the position preference of IR, a unique phenomenon that deserves further molecular characterization.
We found that depletion of the histone methyltransferase SDG725 in rice leads to position-specific alteration of intron retention (IR), a phenomenon that has not been demonstrated previously in any species. Further analyses support the model that H3K36me2/me3 but not H3K36me1 contribute to the IR alteration. As IR plays an important role in regulated gene expression of both plants and animals, the position-specific IR regulation revealed by this study further extends our knowledge regarding the complexity of gene regulatory networks.
Plants of Oryza sativa L. cv. Nipponbare were used in the present study. Seedlings used for RNA extraction and ChIP experiments were grown in artificial growth chambers under a long day (LD) photoperiod (14 h 30 °C: 10 h 28 °C, light: dark). For a rational association between ChIP-seq data and the RNA-seq data, the same set of plant samples was divided into two parts, one part for ChIP-seq library construction, and the other part for total RNA extraction followed by RNA-seq library preparation.
Determination of H3K36 methylation level by mass spectrometry
The 14-day-old rice shoots were subjected to histone preparation based on an extraction method described previously . Then, the prepared histones were separated according to molecular weight with sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) using a standard procedure. The gels were excised at a certain molecular weight and derivatized by In-gel NHS and tested by mass spectrometry (MS) based on a method published recently .
ChIP-seq and RNA-seq library construction and sequencing
We prepared ChIP-seq libraries using 14-day-old rice shoots with the published ChIP protocol . Anti-monomethyl-H3K36 (ab9048), anti-dimethyl-H3K36 (ab9049), and anti-trimethyl-H3K36 (ab9050) were purchased from Abcam, Cambridge, UK. The resulting DNA and input control were subjected to ChIP-seq library construction as previously described . For RNA-seq library construction, total RNA was extracted from 14-day-old rice shoots using TRIzol reagent (Invitrogen). PolyA+ RNA was enriched by selection with two rounds of Oligo(dT)25 Dynabeads (Invitrogen). A strand-specific RNA-seq library was constructed according to a dUTP-based protocol  adapted from Parkhomchuk et al. . Both libraries were sequenced on an Illumina HiSeq2000 instrument to generate 50-bp single-end reads (for ChIP-seq) and 101-bp paired-end reads (for RNA-seq). All related ChIP-seq and RNA-seq data sets are summarized in Additional file 2: Table S9.
Bioinformatics analysis of ChIP-seq data
The sequence tags for H3K36me1/2/3 were aligned against the complete reference genome of Nipponbare (japonica) rice (the MSU Rice Genome Annotation Project, Release 7 ) using Burrows-Wheeler Aligner (BWA) v0.6.1 , with the uniquely mapped reads then extracted with SAMtools . To globally visualize the level of each histone modification along and around rice genes, genes longer than 500 bp were split into 300 blocks, with the upstream and downstream 2-kilobase (kb) regions each split into 100 blocks, respectively. The tag coverage of the aligned nonredundant ChIP-seq reads was then calculated for each block, normalized to 10 million reads, and calibrated to the input. A normalized tag intensity matrix was then obtained with custom scripts, and the average tag intensity was used to construct an aggregation plot. The normalized tag intensity of each gene feature, such as exons and introns, was also calculated as needed. For the heatmaps in Figs. 3 and 4, the x-axis denotes the relative position of a gene from 2 kb upstream of the TSS to 2 kb downstream of the TTS, and the y-axis denotes the normalized coverage of ChIP-seq reads calibrated by the input. The gene body (from TSS to TTS) was split into 300 bins, and the 2-kb upstream TSS and downstream TTS regions were split into 100 bins, respectively; thus, a total of 500 bins for each gene was used to calculate the coverage of ChIP-seq reads. To define the pattern shift of H3K36me2 or H3K36me3 in Fig. 4, we require that the front half (150 bins) of the gene body should have at least 30% of bins (45 bins) with a difference value (725Ri-1 minus WT) greater than 0; at the same time, the latter half should have at least 30% of bins with a difference value smaller than 0. The rest are then considered as not having a pattern shift.
Bioinformatics analysis of RNA-seq data
To calculate the expression abundance of transcripts derived from each gene, we used a series of three programs: Bowtie v1.0.0 , TopHat2 v2.09 , and Cufflinks v2.1.1 [54, 55]. Briefly, the adaptor-removed raw RNA-seq reads were first aligned to Nipponbare (japonica) rice annotated transcripts from Release 7 of the MSU Rice Genome Annotation Project  with Bowtie to estimate insert fragment sizes and standard deviations, which were in turn used as parameter values in TopHat2. TopHat2 was then used to align the paired-end reads to the complete reference genome as mentioned above. Quantification of transcripts and genes, normalized for gene length, was performed with Cufflinks, as represented by fragments per kilobase of exon per million fragments uniquely mapped (FPKM). A differential expression analysis was performed using Cuffdiff, a subpackage of Cufflinks. A cutoff of at least twofold change, FPKM greater than 1 in at least one sample in a sample pair, and a p value smaller than 0.01 was used as the threshold to define DEGs.
Determination of changes in intron retention
To determine the IR and minimize the interference from exons, a new set of introns was acquired so that only introns or intron fragments that do not overlap with any exons were used as corresponding actual introns. In a similar way, we get a new set of exons that do not overlap with any other introns. To determine the IR events, the two exons neighboring an intron must be expressed (i.e., each with an FPKM value greater than 1), and there should be at least three reads supporting the IR events. To evaluate the IR level for a given intron, we introduced the intron retention index (IRI) for an individual intron to quantify its IR level. We obtained the IRI in the following way. The IR level (determined by FPKM) for a given intron was divided by the mean value of the expression level of its neighboring exons (determined by FPKM). To evaluate the IR changes between 725Ri-1 mutants and WT plants, the IRI ratio between the two samples was used, and those exceeding a twofold change were considered as differentially changed (up-regulated or down-regulated) IR events. To compare the location distribution difference between mutant and WT rice, we applied the two-sample Kolmogorov-Smirnov test, which is one of the most useful and general nonparametric statistical methods for comparing two samples without replication, as it is sensitive to differences in both location and shape of the empirical cumulative distribution functions of the two samples .
qRT-PCR for validation of intron retention
Total RNA was extracted from 14-day-old rice shoots using TRIzol reagent (Invitrogen). Reverse transcription was performed using Improm-II Reverse Transcriptase (Promega, Madison, WI, USA) according to the standard protocol provided. Quantitative PCR was performed using the gene-specific primers listed in Additional file 2: Tables S3 and S6.
ChIP-PCR for validation
The 14-day-old rice shoots were used in ChIP-PCR assays as previously described . Anti-dimethyl-H3K36 (ab9049) was purchased from Abcam. ChIP assays were performed according to a previous method . Quantitative PCR was performed to determine the enrichment of immunoprecipitated DNA using a kit from Takara, Otsu, Japan. The gene-specific primers are listed in Additional file 2: Table S6.
Nuclear/cytoplasmic fractionation extraction and quality assay
To determine whether intron-retained transcripts accumulate in the nucleus or cytoplasm, nuclear/cytoplasmic fractionation was performed according to a published method . As quality controls for the fractionation, we evaluated the relative expression abundance of a cytoplasmic marker (Actin1) and nuclear markers (Pri-miR156d, Pri-miR156h, and Pri-miR156b) in the isolated fractions [33, 34]. Strand-specific RNA-seq libraries were constructed using both nuclear and cytoplasmic RNA according to a dUTP-based protocol  adapted from Parkhomchuk et al. .
Cycloheximide treatment and prediction of potential targets of nonsense-mediated mRNA decay
To evaluate the extent that NMD may be involved in the metabolism of intron-retained transcripts, we treated both 725Ri-1 and WT 2-week-old rice plants with cycloheximide (CHX) (which can block the translation process and rescue the NMD-targeted transcripts that would otherwise be degraded ) according to a published method with minor modifications [58, 59]. The same batch of dimethylsulfoxide (DMSO) treatments was used as the control. Then the treated rice samples were used for RNA extraction, followed by RNA-seq library preparation and sequencing as described above. After quality control, the raw sequenced reads were aligned to the rice reference genome as described above, then StringTie was used to assemble transcripts for each sample and combine them together with Cuffcompare, a subpackage of Cufflinks . Then the annotated unique start codon located on the assembled transcript was used for the open reading frame (ORF) prediction. A stop codon was defined as a premature termination codon (PTC) if it was located more than 50 nucleotides upstream of the last exon-exon junction, which is a well-known feature of NMD targets . Lastly, a PTC-containing transcript was defined as a potential NMD target if its accumulative abundance increased more than twofold in 725Ri-1 rice compared to WT plants, with a p value smaller than 0.05, determined by Cuffdiff .
Analysis of splice site strength
MaxEntScan  was used to calculate maximum entropy scores for 9 bp spanning the 5′ (donor) splice sites (3 bp in the exon and 6 bp in the intron) and 23 bp spanning the 3′ (acceptor) splice sites (20 bp in the intron and 3 bp in the exon), respectively.
Wong JJ, Ritchie W, Ebner OA, Selbach M, Wong JW, Huang Y, Gao D, Pinello N, Gonzalez M, Baidya K, et al. Orchestrated intron retention regulates normal granulocyte differentiation. Cell. 2013;154(3):583–95.
Braunschweig U, Barbosa-Morais NL, Pan Q, Nachman EN, Alipanahi B, Gonatopoulos-Pournatzis T, Frey B, Irimia M, Blencowe BJ. Widespread intron retention in mammals functionally tunes transcriptomes. Genome Res. 2014;24(11):1774–86.
Ni T, Yang W, Han M, Zhang Y, Shen T, Nie H, Zhou Z, Dai Y, Yang Y, Liu P, et al. Global intron retention mediated gene regulation during CD4+ T cell activation. Nucleic Acids Res. 2016;44:6817–29.
Guo R, Zheng L, Park JW, Lv R, Chen H, Jiao F, Xu W, Mu S, Wen H, Qiu J, et al. BS69/ZMYND11 reads and connects histone H3.3 lysine 36 trimethylation-decorated chromatin to regulated pre-mRNA processing. Mol Cell. 2014;56(2):298–310.
Shahzad K, Rauf M, Ahmed M, Malik ZA, Habib I, Ahmed Z, Mahmood K, Ali R, Masmoudi K, Lemtiri-Chlieh F, et al. Functional characterisation of an intron retaining K(+) transporter of barley reveals intron-mediated alternate splicing. Plant Biol. 2015;17(4):840–51.
Remy E, Cabrito TR, Batista RA, Hussein MA, Teixeira MC, Athanasiadis A, Sa-Correia I, Duque P. Intron retention in the 5'UTR of the novel ZIF2 transporter enhances translation to promote zinc tolerance in Arabidopsis. PLoS Genet. 2014;10(5):e1004375.
Barbazuk WB, Fu Y, McGinnis KM. Genome-wide analyses of alternative splicing in plants: opportunities and challenges. Genome Res. 2008;18(9):1381–92.
Ner-Gaon H, Halachmi R, Savaldi-Goldstein S, Rubin E, Ophir R, Fluhr R. Intron retention is a major phenomenon in alternative splicing in Arabidopsis. Plant J. 2004;39(6):877–85.
Min XJ, Powell B, Braessler J, Meinken J, Yu F, Sablok G. Genome-wide cataloging and analysis of alternatively spliced genes in cereal crops. BMC Genomics. 2015;16:721.
Wang BB, Brendel V. Genomewide comparative analysis of alternative splicing in plants. Proc Natl Acad Sci U S A. 2006;103(18):7175–80.
Panahi B, Mohammadi SA, Ebrahimi Khaksefidi R, Fallah Mehrabadi J, Ebrahimie E. Genome-wide analysis of alternative splicing events in Hordeum vulgare: highlighting retention of intron-based splicing and its possible function through network analysis. FEBS Lett. 2015;589(23):3564–75.
Li Y, Dai C, Hu C, Liu Z, Kang C. Global identification of alternative splicing via comparative analysis of SMRT- and Illumina-based RNA-seq in strawberry. Plant J. 2017;90(1):164–76.
Gohring J, Jacak J, Barta A. Imaging of endogenous messenger RNA splice variants in living cells reveals nuclear retention of transcripts inaccessible to nonsense-mediated decay in Arabidopsis. Plant Cell. 2014;26(2):754–64.
Boothby TC, Zipper RS, van der Weele CM, Wolniak SM. Removal of retained introns regulates translation in the rapidly developing gametophyte of Marsilea vestita. Dev Cell. 2013;24(5):517–29.
Luco RF, Allo M, Schor IE, Kornblihtt AR, Misteli T. Epigenetics in alternative pre-mRNA splicing. Cell. 2011;144(1):16–26.
Luco RF, Pan Q, Tominaga K, Blencowe BJ, Pereira-Smith OM, Misteli T. Regulation of alternative splicing by histone modifications. Science. 2010;327(5968):996–1000.
Kolasinska-Zwierz P, Down T, Latorre I, Liu T, Liu XS, Ahringer J. Differential chromatin marking of introns and expressed exons by H3K36me3. Nat Genet. 2009;41(3):376–81.
Pajoro A, Severing E, Angenent GC, Immink RGH. Histone H3 lysine 36 methylation affects temperature-induced alternative splicing and flowering in plants. Genome Biol. 2017;18(1):102.
Sui P, Jin J, Ye S, Mu C, Gao J, Feng H, Shen WH, Yu Y, Dong A. H3K36 methylation is critical for brassinosteroid-regulated plant growth and development in rice. Plant J. 2012;70(2):340–7.
Sui PF, Shi JL, Gao XY, Shen WH, Dong AW. H3K36 methylation is involved in promoting rice flowering. Mol Plant. 2013;6(3):975–7.
Liu B, Wei G, Shi J, Jin J, Shen T, Ni T, Shen WH, Yu Y, Dong A. SET DOMAIN GROUP 708. A histone H3 lysine 36-specific methyltransferase, controls flowering time in rice (Oryza sativa). New Phytol 2015. https://doi.org/10.1111/nph.13768.
Kahles A, Ong CS, Zhong Y, Ratsch G. SplAdder: identification, quantification and testing of alternative splicing events from RNA-Seq data. Bioinformatics. 2016;32(12):1840–7.
Chamala S, Feng G, Chavarro C, Barbazuk WB. Genome-wide identification of evolutionarily conserved alternative splicing events in flowering plants. Front Bioeng Biotechnol. 2015;3:33.
Yap K, Lim ZQ, Khandelia P, Friedman B, Makeyev EV. Coordinated regulation of neuronal mRNA steady-state levels through developmentally controlled intron retention. Genes Dev. 2012;26(11):1209–23.
Kalyna M, Simpson CG, Syed NH, Lewandowska D, Marquez Y, Kusenda B, Marshall J, Fuller J, Cardle L, McNicol J, et al. Alternative splicing and nonsense-mediated decay modulate expression of important regulatory genes in Arabidopsis. Nucleic Acids Res. 2012;40(6):2454–69.
Marquez Y, Brown JW, Simpson C, Barta A, Kalyna M. Transcriptome survey reveals increased complexity of the alternative splicing landscape in Arabidopsis. Genome Res. 2012;22(6):1184–95.
Marquez Y, Hopfler M, Ayatollahi Z, Barta A, Kalyna M. Unmasking alternative splicing inside protein-coding exons defines exitrons and their role in proteome plasticity. Genome Res. 2015;25(7):995–1007.
Napoli N, Ghelli R, Brunetti P, De Paolis A, Cecchetti V, Tsuge T, Serino G, Matsui M, Mele G, Rinaldi G et al. A newly identified flower-specific splice variant of AUXIN RESPONSE FACTOR8 regulates stamen elongation and endothecium lignification in Arabidopsis. Plant Cell. 2018. https://doi.org/10.1105/tpc.17.00840.
Jacob AG, Smith CWJ. Intron retention as a component of regulated gene expression programs. Hum Genet. 2017;136(9):1043–57.
Tahmasebi S, Jafarnejad SM, Tam IS, Gonatopoulos-Pournatzis T, Matta-Camacho E, Tsukumo Y, Yanagiya A, Li W, Atlasi Y, Caron M, et al. Control of embryonic stem cell self-renewal and differentiation via coordinated alternative splicing and translation of YY2. Proc Natl Acad Sci U S A. 2016;113(44):12360–7.
Thiele A, Nagamine Y, Hauschildt S, Clevers H. AU-rich elements and alternative splicing in the beta-catenin 3'UTR can influence the human beta-catenin mRNA stability. Exp Cell Res. 2006;312(12):2367–78.
Sun S, Zhang Z, Sinha R, Karni R, Krainer AR. SF2/ASF autoregulation involves multiple layers of post-transcriptional and translational control. Nat Struct Mol Biol. 2010;17(3):306–12.
Wang W, Ye R, Xin Y, Fang X, Li C, Shi H, Zhou X, Qi Y. An importin beta protein negatively regulates MicroRNA activity in Arabidopsis. Plant Cell. 2011;23(10):3565–76.
Ye R, Wang W, Iki T, Liu C, Wu Y, Ishikawa M, Zhou X, Qi Y. Cytoplasmic assembly and selective nuclear import of Arabidopsis Argonaute4/siRNA complexes. Mol Cell. 2012;46(6):859–70.
Hon G, Wang W, Ren B. Discovery and annotation of functional chromatin signatures in the human genome. PLoS Comput Biol. 2009;5(11):e1000566.
Tilgner H, Knowles DG, Johnson R, Davis CA, Chakrabortty S, Djebali S, Curado J, Snyder M, Gingeras TR, Guigo R. Deep sequencing of subcellular RNA fractions shows splicing to be predominantly co-transcriptional in the human genome but inefficient for lncRNAs. Genome Res. 2012;22(9):1616–25.
Yeo G, Burge CB. Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals. J Comput Biol. 2004;11(2–3):377–94.
Liu B, Wei G, Shi J, Jin J, Shen T, Ni T, Shen WH, Yu Y, Dong A. SET DOMAIN GROUP 708, a histone H3 lysine 36-specific methyltransferase, controls flowering time in rice (Oryza sativa). New Phytol. 2016;210(2):577–88.
Hoppmann V, Thorstensen T, Kristiansen PE, Veiseth SV, Rahman MA, Finne K, Aalen RB, Aasland R. The CW domain, a new histone recognition module in chromatin proteins. EMBO J. 2011;30(10):1939–52.
Barski A, Cuddapah S, Cui K, Roh TY, Schones DE, Wang Z, Wei G, Chepelev I, Zhao K. High-resolution profiling of histone methylations in the human genome. Cell. 2007;129(4):823–37.
Wang Z, Zang C, Rosenfeld JA, Schones DE, Barski A, Cuddapah S, Cui K, Roh TY, Peng W, Zhang MQ, et al. Combinatorial patterns of histone acetylations and methylations in the human genome. Nat Genet. 2008;40(7):897–903.
Kuntimaddi A, Achille NJ, Thorpe J, Lokken AA, Singh R, Hemenway CS, Adli M, Zeleznik-Le NJ, Bushweller JH. Degree of recruitment of DOT1L to MLL-AF9 defines level of H3K79 di- and tri-methylation on target genes and transformation potential. Cell Rep. 2015;11(5):808–20.
Li Y, Mukherjee I, Thum KE, Tanurdzic M, Katari MS, Obertello M, Edwards MB, McCombie WR, Martienssen RA, Coruzzi GM. The histone methyltransferase SDG8 mediates the epigenetic modification of light and carbon responsive genes in plants. Genome Biol. 2015;16:79.
Yu Y, Dong A, Shen WH. Molecular characterization of the tobacco SET domain protein NtSET1 unravels its role in histone methylation, chromatin binding, and segregation. Plant J. 2004;40(5):699–711.
Chen J, Gao J, Peng M, Wang Y, Yu Y, Yang P, Jin H. In-gel NHS-propionate derivatization for histone post-translational modifications analysis in Arabidopsis thaliana. Anal Chim Acta. 2015;886:107–13.
Ding B, Zhu Y, Bu ZY, Shen WH, Yu Y, Dong AW. SDG714 regulates specific gene expression and consequently affects plant growth via H3K9 dimethylation. J Integr Plant Biol. 2010;52(4):420-30. https://doi.org/10.1111/j.1744-7909.2010.00927.x.
Zhong S, Joung JG, Zheng Y, Chen YR, Liu B, Shao Y, Xiang JZ, Fei Z, Giovannoni JJ. High-throughput Illumina strand-specific RNA sequencing library preparation. Cold Spring Harb Protoc. 2011;2011(8):940–9.
Parkhomchuk D, Borodina T, Amstislavskiy V, Banaru M, Hallen L, Krobitsch S, Lehrach H, Soldatov A. Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Res. 2009;37(18):e123.
Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S, et al. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013;6(1):4.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, Genome Project Data Processing S. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.
Roberts A, Trapnell C, Donaghey J, Rinn JL, Pachter L. Improving RNA-Seq expression estimates by correcting for fragment bias. Genome Biol. 2011;12(3):R22.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.
Arnold TB, Emerson JW. Nonparametric goodness-of-fit tests for discrete null istributions. R J. 2011;3(2):34–9.
Carter MS, Doskow J, Morris P, Li S, Nhim RP, Sandstedt S, Wilkinson MF. A regulatory mechanism that detects premature nonsense codons in T-cell receptor transcripts in vivo is reversed by protein synthesis inhibitors in vitro. J Biol Chem. 1995;270(48):28995–9003.
Jeong HJ, Kim YJ, Kim SH, Kim YH, Lee IJ, Kim YK, Shin JS. Nonsense-mediated mRNA decay factors, UPF1 and UPF3, contribute to plant defense. Plant Cell Physiol. 2011;52(12):2147–56.
Shen Y, Wu X, Liu D, Song S, Liu D, Wang H. Cold-dependent alternative splicing of a Jumonji C domain-containing gene MtJMJC5 in Medicago truncatula. Biochem Biophys Res Commun. 2016;474(2):271–6.
Maquat LE. Nonsense-mediated mRNA decay. Curr Biol. 2002;12(6):R196–7.
We thank Dr. Wenhui Shen for his critical reading of the manuscript. Thanks also go to Dr. Li Jin and Hillary Sussman for suggestions on the data analysis. We thank Genergy Biotechnology (Shanghai) Co., Ltd. for the deep sequencing service.
This work was supported by the National Key Basic Research Program of China (973 program 2015CB943000 to T.N.) and the National Natural Science Foundation of China (31370752, 31570315, and 91519308 to A.D.; 31471192, 31521003, and 31771336 to T.N.).
Availability of data and materials
The RNA-seq and ChIP-seq raw data are summarized in Additional file 2: Table S8 and can be found at the National Center for Biotechnology Information (NCBI) Sequence Read Archive under accession number SRP136689.
Ethics approval and consent to participate
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. The quantification of H3K36 methylation in 725Ri-1. and WT plants. Figure S2. Differential gene expression between 725Ri-1. and WT rice. Figure S3. Accumulative plot of up-regulated (725Ri-1._IRI_up) and down-regulated (725Ri-1._IRI_down) IR events with different intron length coverage. Figure S4. Distribution of up-regulated (red.) and down-regulated (blue.) intron retention (IR) events between 725Ri-1. and wild-type (WT) rice (a) and between 708Ri-1. and WT (b). Figure S5. Schematic diagram illustrating the way to transform the H3K36 methylation changes between the 725Ri-1. and wild-type plants at individual gene level. Figure S6. Box plot for levels of H3K36 methylations at IRI-up or IRI-down introns. Figure S7. qPCR validation of quality of cellular fractionation in wild-type rice. Figure S8. Box plot for expression level of transcripts with premature termination codon. Figure S9. Box plot of gene expression before and after cycloheximide (CHX) treatment in either wild-type or 725Ri-1. rice. Figure S10. H3K36 methylations associate with gene expression levels. Figure S11. Distribution of up-regulated (red.) and down-regulated (blue.) intron retention (IR) events in genes with no obvious changes (≤ 2-fold) in total levels of H3K36me2 (a) or H3K36me3 (b). Figure S12. Box plot for intron length in two groups of introns by distinct methods in calculating the degree of intron retention. Figure S13. Box plot for GC percentage in two groups of introns. Figure S14. Box plot for maximum entropy score of 5′ splice site (a) and 3′ splice site (b) in two groups of introns. Figure S15. Comparison of H3K36me2/me3 distribution across gene body between rice (a) and Arabidopsis (b). (DOCX 2596 kb)
Table S1. Summary of the RNA-seq libraries. Table S2. Alternative splicing events discovered between rice RNA-seq data sets in this study compared with annotation by using SplAdder. Table S3. Primers used for validating intron-retention events. Table S4. The numbers of IRI changed introns and genes. Table S5. Summary of the ChIP-seq libraries. Table S6. Primers used to validate intron retention in the same gene by qRT-PCR and those used to validate H3K36me2 by ChIP-PCR. Table S7. Number of retained introns for previously identified IRI changed introns in nucleus or cytoplasm. Table S8. Expression changes of annotated rice splicing factors in our data. Table S9. Detailed information for ChIP-seq and RNA-seq data sets. (DOCX 43 kb)
About this article
Cite this article
Wei, G., Liu, K., Shen, T. et al. Position-specific intron retention is mediated by the histone methyltransferase SDG725. BMC Biol 16, 44 (2018). https://doi.org/10.1186/s12915-018-0513-8
- Intron retention
- Histone H3K36 methylation