SF3B1 haploinsufficiency leads to formation of ring sideroblasts in myelodysplastic syndromes

Whole exome sequencing technology has been instrumental in the identification of somatic mutations involving genes important in RNA splicing.^1–4  The first and most frequently identified mutation in hematologic malignancies involves the SF3B1 gene and was found in 60%-80% of refractory anemia with ring sideroblasts (RARS) and RARS associated with thrombocytosis (RARS-T) patients. Both RARS and RARS-T are subtypes of myelodysplastic syndromes (MDSs) and MDS/myeloproliferative neoplasm (MDS/MPN) overlap diseases. SF3B1 is a crucial component of the spliceosomal U2snRNP complex, which participates in normal RNA splicing.⁵  Intron splicing that culminates in the production of the mature mRNA product is catalyzed by spliceosomes. The major spliceosomes consist of U1, U2, U4/U6, and U5snRNPs and carry out splicing of major class or U2-type introns. Minor spliceosomes, which consist of U11, U12snRNPs, U4atac/U6atac and U5 snRNPs, carry out the splicing of minor class introns or U12-type introns. RNA splicing is crucial in maintaining normal genetic diversity.⁶  Because the fidelity of the final protein product depends highly on an intact and accurate RNA splicing machinery, somatic genetic errors in key components of the splicing pathway can potentially lead to the formation of dysfunctional proteins that can ultimately predispose to diseases. Indeed, emerging discoveries of genetic polymorphisms or mutations altering regulatory sequences or producing splice variants have been associated with certain types of hereditary diseases and neoplastic disorders.⁷ 

The presence of 15% or more ring sideroblasts (RSs) of nucleated erythroid precursors in the BM is a cardinal feature of RARS and RARS-T, even though RSs are also seen in 25% of patients with other subtypes of MDS, but rarely in other hematologic disorders.^8–10  RS occurs because of abnormal localization of crystalline iron deposits in the mitochondria of erythoid precursors. Even though there is evidence of dysfunction in mitochondrial metabolism, the pathogenesis of RS remains unclear.^11,12  The discovery of SF3B1 mutations in RARS and RARS-T has opened a new area of investigation in these diseases, and this article focuses on experimental evidence linking dysfunction in SF3B1 and formation of RS in RARS/RARS-T.

BM and peripheral blood cells were obtained with written informed consent in accordance with the Declaration of Helsinki and approved by the Institutional Review Board of the Cleveland Clinic. Diagnosis was assigned according to 2008 World Health Organization classification criteria. We studied a total of 456 patients. Overall, the MDS group includes: RARS, n = 19; refractory anemia with excess blasts 1/2, n = 34; refractory cytopenia with unilineage dysplasia/refractory cytopenia with multilineage dysplasia, n = 23; and MDS-U/5q, n = 17. The MDS/MPN group includes: RARS-T, n = 26; chronic myelomonocytic leukemia 1/2, n = 59; and MDS/MPN-U, n = 20. The acute myeloid leukemia (AML) group includes: secondary AML, n = 51; and primary AML, n = 43. Others includes: MPN, n = 69; mast cell disease, n = 31; and paroxysmal nocturnal hemoglobinuria/aplastic anemia/T-cell large granular lymphocytic leukemia, n = 64. The median age at sample collection was 52 years (range, 13-77 years).

Genomic sequencing was performed on coding regions for SF3B1 (exons 13-16) SF3B14 (all exons), SF3B4 (all exons), and DYRK1A (all exons). Primer sequences and conditions used can be supplied on request. Bidirectional sequencing was performed by standard techniques using an ABI 3730xl DNA analyzer (Applied Biosystems). Somatic nature of novel mutations was confirmed using germline DNA from CD3⁺ cells. All mutations were scored as pathogenic on the basis of the observation that they were not detected in normal samples and were not found in published SNP databases (dbSNP, http://www.ncbi.nlm.nih.gov/projects/SNP) and/or they were not reported as SNPs in previous publications.

One million BM mononuclear cells were first subjected to fixation in 3.75% glutaraldehyde, 0.1M sodium cacodylate, 6% sucrose (pH 7.2-7.4) for a period of 4 hours to overnight. Specimens were post-fixed in 1% osmium tetroxide and 0.1M cacodylate (pH 7.2-7.4) for 1 hour. After each fixation step, specimens were washed twice in cacodylate buffered sucrose, 10 minutes each. After the final wash, specimens were dehydrated in ascending graded alcohols, followed by immersion in propylene oxide, then passed through ascending graded mixtures of propylene oxide and EMBed 812 epoxy resin, and finally embedded in pure epoxy resin. Specimen blocks were polymerized overnight at 60°C. Specimens were initially evaluated by light microscopy of 1-μm plastic sections stained in a combination of toluidine blue and basic fuchsin. Thin sections from selected tissue blocks were cut at 60-80 nm and supported on 200 mesh copper grids. Sections were stained in uranyl acetate and lead citrate and examined using a Philips CM-12 electron microscope. Electron micrographs were taken using a Gatan Orius CCD digital camera.

Virus containing supernatant of 293T cells was used to infect 5 × 10⁵ K562 cells. Five different shRNAs were constructed by Sigma-Aldrich. To assess the effect of shRNAs, cells were harvested 2 days after puromycin selection and SF3B1 mRNA level was measured by RT-PCR. Total RNA was isolated using the High Pure RNA Isolation kit (Roche Diagnostics). For quantitative RT-PCR, 0.5 μL of the complementary DNA was added to SYBR@advantage qPCR mix (Clontech Laboratories) and amplified in a Light Cycler 480 instrument (Roche Diagnostics). The U2- and U12-dependent introns were selected from the online database (http://genome.crg.es/datasets). The RT-PCR primers spanning the exon-intron junctions of several genes were designed using PrimerQuest (www.idtdna.com). And their specificity was assessed in silico (http://genome.ucsc.edu) and by standard PCR. A total of 40-50 PCR cycles were performed in a 2-step cycling procedure with an initial denaturation at 94°C for 3 minutes, subsequently at 94°C for 15 seconds, and 60°C for 30 seconds. Data were analyzed based on relative δ C_t values using GAPDH as control gene. Results are shown as mean ± SD of duplicate experiments.

A total of 100 000 BM cells derived from 4 healthy persons were plated in duplicate in 1 mL methylcellulose with cytokines (StemCell Technologies) in 35-mm culture plates for 2 weeks at 37°C with 5% CO₂. Meayamycin was added at the beginning of the culture at a range of concentrations (1, 2, 5, 10, and 50nM). Colony-forming cells were harvested, washed 5 times in PBS, and spotted on cytospin slides before Prussian blue staining.

Meayamycin was synthesized according to the literature¹³  and purified on a C18 reverse-phase preparative HPLC. The fractions containing meayamycin were concentrated under vacuum to afford 16.3 mg of solid meayamycin. Using this solid, a 10mM stock solution of meayamycin in DMSO was prepared. A 10μM solution of meayamycin in DMSO was prepared through serial dilutions and stored at a −20°C.

BM aspirates from Sf3b1 heterozygous¹⁴  and wild-type (WT) C57BL/6 mice were kindly provided by Dr Koseki (RIKEN) and were subjected to Prussian blue staining using standard pathology procedures for the detection of RS.

RNA sequencing (RNA-Seq) was performed by Otogenetics Corporation (Tucker); 100-bp paired-end RNA-Seq reads were mapped to the hg19 RefSeq human transcriptome and spliceosome by DNAnexus (http://dnanexus.com) using a Bayesian method,¹⁵  in which a read was mapped when its posterior probability of mapping exceeded 0.9. These filtered posterior probabilities were summed to generate fractional read counts per gene and per exon, with probabilities from splice-junction spanning reads counted for each relevant exon. We used rounded gene and exon read counts as inputs for our differential expression analyses.

We used the R package DESeq¹⁶  version 1.6.1 (http://www.bioconductor.org/packages/release/bioc/html/DESeq.html) and R version 2.14.1, to perform differential gene expression analysis of mutant versus normal samples. DESeq uses a negative binomial (NB) distribution to model the gene read counts and shrinkage estimators to estimate the per-gene NB dispersion (square of the coefficient of biologic variation) parameters. Specifically, we estimated the per-gene dispersion parameters using a shrinkage approach with the default maximum sharing mode and default parametric mean-dispersion relationship, where raw dispersion values were calculated using the samples in the mutant group. Before testing, we dropped all genes with read counts below the 40th percentile to improve testing power while maintaining type I error rates. We used the nbinomTest function, which does conditional testing assuming gene counts follow an NB distribution with fitted dispersions, to estimate P values and adjusted P values obtained using the Benjamini-Hochberg method¹⁷  for each filtered gene under the null hypothesis of common expression intensity across groups. Genes with adjusted P < .05 were declared significant. Moderated logarithm base 2-fold changes (mutant/normal) were estimated using a variance stabilizing transformation as implemented in the getVarianceStabilizedData function of DESeq.

We used the R package DEXSeq Version 1.1.4 (http://www.bioconductor.org/packages/release/bioc/html/DEXSeq.html), to perform differential exon usage analysis of mutant versus normal samples. DEXSeq uses an NB distribution to model the exon read counts and shrinkage estimators to estimate the per-exon NB dispersion parameters. Specifically, we estimated the per-exon dispersion parameters using a shrinkage approach with the maximum sharing mode and parametric mean-dispersion relationship, where raw dispersion values were calculated using the samples in the mutant group. We defined a testable exon as one that had a total sum of at least 10 mapped reads across samples and was in a gene with no more than 70 exons. Before exon usage testing, we dropped any exons that were not testable or were in genes with less than 2 testable exons to improve testing power while maintaining type I error rates. We used the testGeneForDEU function, which compares deviances from generalized linear model fits (assuming NB likelihood) to a χ² reference distribution, to estimate P values and adjusted P values obtained using the Benjamini-Hochberg method for each exon under the null hypothesis of common usage across groups. Exons with adjusted P < .05 were declared significant. Logarithm base 2-fold changes (mutant/normal) for each exon were estimated using the function estimatelog2FoldChanges.

We used the R package goseq¹⁸  Version 1.6.0 (http://www.bioconductor.org/packages/release/bioc/html/goseq.html) to perform over- and under-representation analysis (ORA) of gene sets taking length bias into account. We conducted separate ORAs for the gene-level and exon-level differential expression (DE) findings. ORA analysis requires a definition of significance for each gene member of the gene set. We defined this similarly for ORA of the gene-level and exon-level DE findings: (1) a gene with an adjusted P value of < .5 from the gene-level DE analysis and (2) a gene with at least one exon with an adjusted P value from the exon-level DE analysis < .5, respectively. Goseq requires a probability weighting function, which is used to adjust the gene set results to reflect the lengths of genes in the gene sets, reducing the common bias in gene set analyses because of longer genes having higher read counts and subsequently a higher chance of differential detection. For ORA of the gene-level DE findings, we defined length as the sum of the lengths of the observed exons within each gene, where the exon lengths were obtained from DNAnexus. For ORA of the exon-level DE findings, we defined length as the median exon length of detected exons within the specified gene, where exon lengths were obtained from DNAnexus. We used the MSigDB database¹⁹  Version 3.0 (http://www.broadinstitute.org/gsea/msigdb/index.jsp), gene set collections c2 through c5 in our ORA. For each MSigDB collection of gene sets, we ran the goseq function with the Wallenius approximation to estimate P values for over-representation and under-representation of significant genes (where significant is defined earlier in this paragraph) for each gene set within the collection. Adjusted P values were calculated on the gene-set P values per collection using the Benjamini-Hochberg method to control for the number of gene sets tested within a collection. These adjustments were made separately for the P values from the over-represented tests and for the P values from the under-represented tests. Any gene sets with adjusted P values < .05 from either the over- or under-represented tests were declared significant.

Overall survival was measured from the day of initial sampling to death from any cause (patients lost to follow-up were censored) or last follow-up and was summarized using the Kaplan-Meier method. Results were analyzed for data collected as of December 2011. All P values were 2-sided, and P values ≤ .05 indicated statistical significance. Statistical analyses were performed using JMP8 (SAS Inc). RNA-Seq data statistical analysis is discussed extensively in “RNA-Seq analysis.”

We first examined SF3B1 mutational status in our cohort of MDS and other hematologic malignancies (n = 456) using Sanger sequencing. Consistent with earlier studies,^1–3  we found 50 cases with somatic mutations in SF3B1; 68% of RARS and 81% of patients with RARS-T harbored mutations, whereas SF3B1 mutations were infrequently found in higher-risk MDS and other hematologic conditions (supplemental Figure 1A-D, available on the Blood Web site; see the Supplemental Materials link at the top of the online article). The K700E was the most frequent mutation detected in SF3B1 among patients with RARS and RARS-T (34%). Further, SF3B1 did not predict for progression as underlined by the low prevalence of mutations in patients with secondary AML (5.9%) and de novo AML (4.7%; supplemental Figure 1E). Analysis of overall survival comparing mutant versus WT patients showed longer survival in patients with SF3B1 mutations (supplemental Figure 3). Interestingly, we next found 2 non-MDS cases in which SF3B1 mutations were detected: one with paroxysmal nocturnal hemoglobinuria (K666Q)²⁰  and one with post-polycythemic myelofibrosis (K700E).¹  Both displayed RS (17% and rare, respectively). This finding led us to investigate more closely the association of SF3B1 mutations and RS. In conclusion, SF3B1 is mainly associated with RARS and RARS-T MDS and MDS/MPN subtypes. In a total of 45 patients, 11 typical RARS and RARS-T patients with more than 15% RS were WTs for SF3B1. Thus, we sought to investigate whether SF3B1 might accompany distinctive features.

Given the high frequency of SF3B1 mutations in patients with RARS,^1–3,21  we then investigated the association of spliceosomal protein mutations with altered mitochondrial iron distribution manifested as RS. We analyzed ultrastructural differences of BM erythroid precursors between SF3B1 mutant and WT patients with the intent of evaluating whether SF3B1 mutations might alter iron distribution. Transmission electron microscopy was performed on BM cells derived from an RARS patient carrying an SF3B1 mutation (K700E) and an RARS patient WT for SF3B1. Abundant foci of iron deposits were observed around the nucleus of BM cells in the mutant case compared with scanty deposits in the WT patient. Although Prussian blue staining is traditionally used to identify RS in BM samples, it does not distinguish between SF3B1 mutant and WT cases (Figure 1 top right panel). The SF3B1 mutations are associated with more dense iron particles. It is possible that the topographic configuration of the iron might predict the presence of SF3B1 mutations.

We then investigated the role of SF3B1 in pre-mRNA splicing of endogenous human genes. First we knocked down the level of SF3B1 mRNA in cells using RNAi. We applied shRNA lentiviral infection in K562 cells using 5 different shRNAs specific for the human SF3B1 gene. We performed quantitative RT-PCR to measure the level of knockdown. SF3B1 mRNA knockdown of 50% was consistently achieved on one shRNA clone. Using this selected shRNA, we quantified the ratio of spliced and unspliced RNA in 7 and 3 exemplary U2- and U12-introns. All examined U2-type introns were less efficiently spliced, whereas U12-introns were not significantly affected (Figure 2A; Table 1). We then analyzed the phenotype of K562 cells with decreased expression of SF3B1 mRNA. K562 cells were first cocultured with hemin, an inducer of hemoglobinization and erythroid differentiation, and then infected with different shRNAs specific for SF3B1. A 50% reduction of SF3B1 mRNA level was again achieved. Transfectant cells were sent for morphologic analysis and Prussian blue staining. Increased mitotic figures were noted, but RSs were not observed. Transmission electron microscopy of the transfectant cells also showed no iron deposits compared with controls (ie, K562 cells infected with a vector only; data not shown).

To further explore the relationship between SF3B1 mutations and RS, we used meayamycin, a chemical analog of FR901464.^13,22  FR901464 is a potent natural product that inhibits pre-mRNA splicing via SF3b factor binding.^23,24  We performed human colony-forming unit cell assay using BM mononuclear cells from 4 healthy subjects without and with different concentrations of meayamycin (2, 10, and 50nM). The colonies were spotted on cytospin slides and stained with Prussian blue. As shown by the representative pictures, meayamycin consistently induced RS in all 4 healthy subjects compared with controls (Figure 2B). The percentages of RS per each meayamycin concentration are shown in supplemental Figure 2.

To next examine the role of SF3B1 in RS formation, we analyzed the phenotype of a heterozygous mouse model for Sf3b1 (Sf3b1^+/−) insufficiency in vivo.¹⁴ Sf3b1^+/− mice showed only minimal differences in blood counts, which includes a mild lymphopenia compared with WT C57BL/6 mice (21.9 × 10²/μL ± 2.9 vs 42.3 × 10²/μL ± 11.6; n = 3). Indeed, hemoglobin levels (8.5 g/dL ± 0.2 vs 7.9 g/dL ± 0.2) and platelets counts (60.2 × 10⁴/μL ± 13.4 vs 59.8 × 10⁴/μL ± 5.4) were comparable between heterozygous and WT mice. C57BL/6 mice served as controls because Sf3b1^+/− mice were backcrossed to this mouse strain. BM aspirates were provided by Dr Koseki (RIKEN) and stained for Prussian blue. We found numerous RS in Sf3b1^+/− compared with C57BL/6 mice as shown in Figure 2C (right panel).

To further understand the functional consequences of SF3B1 mutations, we used RNA-Seq to investigate the differential exon usage and gene expression levels in BM cells of 2 SF3B1 mutant RARS patients. Using a false discovery rate level of 0.05, we tested 81 564 exons in 9069 genes and found that 423 exons derived from 350 genes were differentially used in mutants, compared with a representative healthy donor. A total of 67% of these exons had higher usage in mutant versus normal, and 56 genes were significantly affected in more than 1 exon (Table 2). Of note, 4 genes relevant to MDS pathophysiology, namely, ASXL1, CBL, EZH1, and RUNX3, were found to be alternatively spliced in at least one exon (Figure 3). Gene level analysis conducted on 12 870 reference genes showed that none of these 4 genes had differences in expression level. However, we identified 130 differentially expressed genes, of which 94% had lower expression in SF3B1 mutant patients (supplemental Table 1). Gene set analysis was also performed on 11 773 genes showing many gene sets related to myeloma or MDS (supplemental Table 2). In addition, a cytokine gene set, cytokine production, was the 16th most significant among the c5 collection of gene sets. This observation supports prior findings that cytokines play an important role in MDS pathogenesis.