The transcriptomic landscape of normal and ineffective erythropoiesis at single-cell resolution

Erythropoiesis is a complex process in which bone marrow progenitors multiply and differentiate into mature, enucleated red blood cells. Humans produce 2.3 × 10⁶ red blood cells each second, and this can increase from fivefold to 10-fold when needed, as documented based on the reticulocyte count. To allow high and adaptable output while maintaining cell integrity, erythropoiesis needs to be tightly regulated, efficient, and effective. Should excessive numbers of erythroid precursors die while maturing in the bone marrow, anemia ensues, and erythropoiesis is termed ineffective. Diamond Blackfan anemia (DBA), thalassemia, and the anemia of myelodysplastic syndromes (MDSs)¹ are examples of ineffective erythropoiesis. However, despite previous studies by our group and others,^2-12 the reasons why erythroid cells die is unclear.

In contrast, the sequential events of effective erythropoiesis are well understood. As burst forming unit–erythroid (BFU-E; earliest erythroid progenitors) mature to colony forming unit–erythroid (CFU-E) and proerythroblasts, their cell surface transferrin receptor (TfR1; CD71) is upregulated, allowing an influx of iron and upregulation of heme synthesis. Heme (a toxic chelate) induces the transcription and translation of globin by removing inhibitors Bach1 and HRI, respectively. At the time when iron is plentiful and heme synthesis is robust yet globin is insufficient, early erythroid precursors depend on feline leukemia virus C receptor^13,14 to export unneeded intracellular heme and avoid heme- and reactive oxygen species (ROS)-mediated damage via apoptosis and ferroptosis.^3,15-17 Later in erythroid differentiation, many cell-intrinsic mechanisms adjust red blood cell size (mean corpuscular volume) in response to heme or globin availability.^18,19 

Erythropoiesis is commonly tracked using flow cytometry because maturing erythroid cells sequentially express cell surface markers CD71, CD36, glycophorin A, and band 3, and these markers correlate well with developmental stage and cell function.^20-22 Recent data, however, suggest that differentiation, including erythropoiesis, is a continuum rather than discrete steps.^23-25 Ostensibly, each cell faces obstacles during maturation, and its survival depends on its internal programming and ability to overcome successive challenges or compensate for them. Because insults and compensations cannot be separately analyzed after flow-cytometric sorting or transcriptional clustering, there is a compelling rationale for single-cell methods.^12,26,27 

Therefore, to delineate the shared mechanisms that lead to the death of maturing erythroid cells, we studied single marrow cells from patients with MDS-5q, patients with DBA, and healthy individuals using CITE-seq (cellular indexing of transcriptomes and epitopes by sequencing; antibody barcoding of cell surface proteins)²⁸ and single-cell RNA sequencing (scRNAseq). We assigned cell stages using traditional flow cytometry surface markers, overlaid transcriptional pseudotime determinations, defined healthy vs unhealthy developmental trajectories, and analyzed gene expression as differentiation proceeded.

In addition, our observations help address an enigma in patients with MDS at low risk (at presentation, often only 50%-75% of cells derive from the mutant clone, yet anemia is prominent). Why residual healthy (unmutated) cells cannot expand and compensate is unknown. We developed methods to independently track precursors with and without chromosome 5q, resolve insult from compensation, and derive new insights into the cell-extrinsic regulation of erythropoiesis.

To investigate erythroid differentiation during normal and ineffective erythropoiesis, lineage-depleted marrow cells were analyzed via flow cytometry and processed for antibody barcoding and scRNAseq either directly (day 0) or after 3 or 6 days of erythroid expansion culture (Figure 1; supplemental Figure 1). Marrow samples were collected after obtaining written informed consent with an institutional review board–approved protocol. We generated a scRNAseq data set from 4 patients with ineffective erythropoiesis and anemia (2 patients with DBA and 2 patients with MDS-5q) and 7 healthy controls (Table 1) and sequenced 74 214 individual marrow cells after excluding doublets.³⁰ Cells containing >15% mitochondrial reads were considered to be dead so were excluded (4869 removed and 69 345 retained). Moreover, genes expressed in <20 cells were filtered to assure robust statistical analysis. We then retained the 7500 most variable genes across all cells (pp.highly_variable_genes function). At least 1 control sample was included in each experiment, and some samples were included in multiple experiments to aid comparisons of independent experiments and to demultiplex samples (supplemental Table 1).

By day 6, all cultures contained ample early erythroid precursors (cluster Ery1; Figure 1D). Expression of some globin genes and TFRC were increased in culture as expected (supplemental Figure 2A-B).^31,32 Importantly, the transcriptional program of cultured cells reflected that of cells in vivo, and principal component analysis showed only minimal differences (supplemental Figure 2C-D), justifying the combined analyses.

RNA velocity calculations of the ratios of spliced and unspliced RNAs were obtained with the pp.moments and tl.velocity functions of scVelo.³³ We generated pseudotime trajectories using the tl.velocity_pseudotime function, which uses both spliced and unspliced transcripts. Multiplexed samples were demultiplexed based on natural genetic variation without assistance of donor genome references^34-36 and then identified by checking for known mutations, attenuated RPS14 gene expression (5q deletion [5q^–]), or sex-specific gene expression. Gene set scores were calculated using the Scanpy³⁷ score_genes function. Default parameters were used in all calculations, unless noted otherwise. Eleven cell-type scores were calculated using panels of marker genes,^38,39 and their expression was evaluated over time in culture (Figure 1B-E; supplemental Figure 3). Scores were also calculated for cell cycle phase using the gene sets and methods of Satija et al.⁴⁰ Cells were clustered via Leiden community detection,⁴¹ using batch-corrected nearest neighbors.⁴² Differentially expressed genes were identified using the Wilcoxon rank-sum test with Benjamini-Hochberg correction (P < .05); then, gene ontology (GO)-term enrichment analysis was performed using the Enrichr method,⁴³ which calculates the odds ratio based on the deviation from expected rank. Experiments were analyzed both separately and combined after normalization and scaling (Scanpy functions), using multiple algorithms to ensure individual samples or algorithms did not bias results. We used uniform manifold approximation and projection⁴⁴ visualization as an atlas for all cell types, diffusion maps for 2-dimensional visualization of erythroid differentiation, and pseudotime scatter plots to visualize changes in individual genes or pathways over differentiation. See supplemental Methods for additional details, including how independent data sets were integrated and analyzed.

We evaluated multiple approaches to order cells, including several pseudotime methods (supplemental Figure 4), and opted for the velocity pseudotime method of scVelo³³ because it incorporates both spliced and unspliced messenger RNA (mRNA) in its calculations. The large abundance of incompletely spliced mRNAs during erythropoiesis^45,46 provided increased quantities of RNA, which resulted in a linear pseudotime that was consistent across different experiments. To show that our pseudotime ordering appropriately reflected sequential differentiation, we constructed individual pseudotime plots at baseline (day 0) and at culture days 3 and 6. Cell progression was as expected (supplemental Figures 1 and 5). Moreover, cells from patients with DBA and those with MDS progressed more slowly through pseudotime than cells from healthy individuals. This reinforced the need to compare accurately staged cells rather than cells at equivalent times in culture. Next, we used the expression and loss of traditional flow cytometry differentiation markers,²⁸ CD117 (cKit), CD36, CD71, and CD235a (GlyA) to further validate the pseudotime order (Figure 2A-D; supplemental Figure 6A-E). Importantly, this also allowed us to assign erythroid stages BFU-E through orthochromatic erythroblast (stages B-E) to the cells along the pseudotime order based on the known stage-specific regulation of these well-studied cell surface proteins.^21,47,48 We then confirmed the assigned stages using other published data, such as the expression patterns of HBB, CD34, TRIB2, and CA2, and the increased frequency of cells in S phase during the CFU-E stage (Figure 2E-H; supplemental Figure 6F).^{21,26,27,49,50} Unexpectedly, the frequency of cells in S, M, and G phase at each stage of differentiation was identical among samples (Figure 2I). Thus, the rate of cell cycling appeared intrinsic to the differentiation stage and not the specific disease. This suggests that the rate of cell division is linked to developmental processes occurring at particular stages of differentiation and that differences observed in population-based cell cycle studies could reflect differences in the number of cells at each stage present in the samples and not the disease pathophysiology. The number of distinct genes expressed during erythropoiesis peaks at the BFU-E stage and then steadily declines during terminal differentiation (Figure 2J).

Although alternative splicing of mRNA is known to be a key regulatory mechanism,^51-53 intron retention is also prevalent throughout erythroid differentiation. Because unspliced mRNA precedes spliced mRNA, the ratio of unspliced to spliced mRNAs within each cell, or RNA velocity,³³ can predict developmental trajectories in many cell types. Interestingly, this ratio increased as erythroid precursors matured to the point that the RNA velocity algorithm predicted a reversed trajectory of differentiation for erythroid cells after the CFU-E stage but not for any other cell lineage (Figure 3A). This finding has previously been reported in fetal erythropoiesis.⁵⁴ Along with the stage-specific changes in intron retention, we noted gene-specific changes in intron retention. Genes such as DDX39B, CA2, KLF1, and TFRC showed peak intron retention earlier in differentiation, whereas intron retention in GATA1, TAL1, and globin gene transcripts (eg, HBB) peaked at the terminal stage of erythropoiesis (Figure 3B). We then analyzed the ratio of unspliced to spliced mRNA in bulk ex vivo and cultured marrow via quantitative polymerase chain reaction and observed similar changes (Figure 3C). Thus, erythroid gene expression has broad posttranscriptional regulation through intron retention in addition to alternative splicing. As messages with retained introns are not likely to be effectively translated, this provides an additional method to regulate protein expression during erythropoiesis.

Erythroid precursor cells from patients with DBA or MDS-5q upregulated stress response, apoptosis-related, and mitochondrial metabolism genes (supplemental Figure 7; supplemental Table 2). ADA was highly upregulated in patient-derived DBA cells, as expected.^55,56 

As a more rigorous and informative approach, we analyzed the stage-specific transcriptome to identify the differences between healthy controls and patients with anemia, which were consistent across all experiments and times. When we visualized the changes in transcriptomes over developmental time, we noted 2 prominent differentiation trajectories that branched at the CFU-E stage and were present in all samples (Figure 4A). The main trajectory (trajectory A) continued through terminal differentiation, whereas the alternative trajectory (trajectory B) ended at or before the proerythroblast stage. Each trajectory included cells from both patients with anemia and healthy individuals (Figure 4A) and were present with many visualizations and algorithms (supplemental Methods), and after down sampling (supplemental Figure 8). Our analysis of published scRNAseq data from a study of marrow CD34⁺ cells isolated from patients with DBA and healthy individuals by Iskander et al¹² also showed these 2 cell-fate trajectories in all samples (supplemental Figure 9). These 2 trajectories were also present but unappreciated in studies of human fetal liver erythropoiesis.⁵⁴ Thus, we concluded that it was a common feature of erythropoiesis.

To understand the differences between cells on trajectory A vs those on trajectory B, we analyzed gene and pathway expression of equivalently staged cells on each trajectory. The top 20 genes upregulated in trajectory-A cells were mostly erythroid differentiation genes and included some cell cycle and mitochondrial oxidative metabolism genes. These genes were poorly expressed in trajectory-B cells (Figure 4B-D; supplemental Table 3). Although many new genes were upregulated when cells progressed from early CFU-E (C1) to late CFU-E (C2) along trajectory A,10-fold fewer new genes were expressed when cells progressed from C1 to C2 along trajectory B (32 vs 348; Figure 4E). Unbiased GO enrichment analysis replicated this observation. Most pathways were downregulated in C2 along trajectory B in healthy donor samples and in samples from patients with DBA or MDS-5q (Figure 4D; supplemental Table 3), indicating failed differentiation.

Of the day-0 (uncultured) B- and C-stage cells from healthy individuals, 24% ± 12% followed trajectory B, whereas twofold or threefold more day-0 cells from patients with DBA or MDS-5q (73% ± 2% and 57% ± 4%, respectively) followed trajectory B (supplemental Table 4). However, culturing allowed more patient cells to follow trajectory A and complete erythroid differentiation (supplemental Table 4). This finding confirms our previous work, showing that culture conditions enable more erythroid precursors from patients with DBA or MDS-5q to survive.³ 

Protein barcoding and pseudotime analysis placed the trajectory branch point during the CFU-E stage, thus, after iron import via CD71 and when heme synthesis is first upregulated (Figures 2 and 5A-B). Despite significant transcriptional differences between cells on trajectory A vs those on trajectory B (Figure 4B-D), the cell surface expression of CD36, CD71, and GlyA on cells as they progressed along both trajectories were sufficiently similar to prevent prospective flow-based resolution of these trajectories (Figure 5B).

Cells on trajectory B showed a rapid decrease in nuclear gene expression and a corresponding increase in the percent of mitochondrial reads (Figure 4F), indicating a loss of membrane integrity. Moreover, these cells had the highest expression of BAX and high levels of nonspecific CITE-seq antibody binding, all factors consistent with trajectory B being a death pathway. A cell’s decision to follow trajectory B (per transcriptomics) occurs well before its ultimate demise (Figure 5A-B). This implies that there is time during early erythroid differentiation when therapeutic interventions could amend cell fate.

Although very few genes were upregulated in trajectory-B compared with those in trajectory-A cells, the highest expressed genes in trajectory-B cells in all samples were ribosomal protein genes. These heme-responsive genes⁵⁹ were significantly upregulated in early CFU-E, the stage at which trajectories A and B initially diverge (Figures 4C and 5C). We also noted a significant increase in heme-responsive gene expression, as defined by Liao et al,⁶¹ in C2 stage trajectory-B cells (Figure 5D), the stage during which heme synthesis intensifies. In addition, pseudotime analysis revealed significant downregulation of GATA1 targets, hallmark heme metabolism, mitotic spindle, and ribosome pathway genes in trajectory-B vs trajectory-A cells (Figure 5C). This is in addition to the downregulated mitochondrial electron transport and mRNA splicing genes revealed by GO analysis (Figure 4D). Because these findings were present in trajectory-B cells in all samples, and because there results are similar to the transcriptional alterations observed in our single-cell analyses of Flvcr1-deleted and Rpl11-haploinsufficient mice,^59,60 the data implicate heme excess in cell demise. The transcriptional alterations in trajectory-B cells are also similar to those seen in in vitro studies of normal human CD36⁺ erythroid cells 15 minutes after heme synthesis is induced by the addition of aminolevulinic acid.⁵⁹ That the transcriptomes of the erythroid cells from healthy individuals on trajectory B resembled the transcriptomes of DBA and MDS-5q erythroid cells on trajectory B was an unexpected finding and suggests that heme regulates normal as well as pathological erythroid maturation.

Next, we cultured healthy marrow cells and cells from a patient with DBA with succinylacetone to inhibit heme synthesis, knowing that this reduces cell death and increases the maturation of red blood cells in DBA samples.³ This resulted in a 24% increase in the frequency of DBA cells on trajectory A (Figure 5E) and resulted in increased expression of GATA1 targets, heme-repressed genes, and BACH1 activity (Figure 5F), as expected from reduced heme production. There was decreased expression of key proapoptotic genes and increased expression of some antiapoptotic genes (Figure 5G), thus implying that erythropoiesis was more effective.

Some DBA and MDS-5q cells can follow trajectory A and survive. To determine how this occurs, we performed GO analyses of erythroid precursor cells at each stage of trajectory A differentiation and compared this with the results of samples from healthy individuals (Figure 6; supplemental Tables 5-6). This allowed us to identify alterations that were common to both DBA and MDS-5q.

Cell cycle, cell cycle mitotic, and dimerization partner, Rb-like, E2F, and MuvB target pathways were significantly dysregulated in DBA cells on trajectory A when compared with healthy cells on trajectory A, indicating that DBA trajectory-A cells were experiencing stress. Ribosomal protein genes (reactome translation pathways) were upregulated in early DBA trajectory-A erythroid cells (Figure 6A; supplemental Table 5), consistent with elevated heme levels.^59,60 p53 signaling pathway, p53 targets, and hypoxia pathway genes were upregulated later. Notably, DBA trajectory-A erythroid cells had significantly reduced expression of heme biosynthesis genes and related pathways (heme biosynthesis, porphyrin and chlorophyll metabolism, and α-hemoglobin stabilizing enzyme pathways) beginning at the BFU-E stage. Although GATA1 was transcriptionally upregulated, GATA1 targets were not uniformly upregulated (heme biosynthetic gene targets were downregulated whereas other targets such as globin genes were upregulated). Thus, surviving DBA cells (trajectory A) downregulate heme production and upregulate globin transcription. This might offset their slowed protein (globin) translation from ribosomal haploinsufficiency, mitigate heme toxicities, and permit survival.

Unbiased GO analysis of erythroid precursor cells on trajectory A from patients with MDS-5q showed similar alterations and only minor differences (Figure 6B; supplemental Table 6), indicating that the shared pathophysiology between DBA and MDS-5q results in shared compensatory responses. Interestingly, only RPS19, but not other ribosomal proteins, was upregulated in C1 trajectory-A MDS-5q cells. The upregulation of HIF1A targets, hypoxia, and DNA damage response genes in trajectory-A DBA and MDS-5q cells would counteract the ROS produced by excess heme and further promote cell survival. These changes were not seen in trajectory-B cells.

The ribosomal haploinsufficiencies present in DBA and MDS-5q are sometimes associated with increased p53 activity,^4,62-64 and the TP53 network and p53 signaling pathway and DNA damage response pathway were commonly upregulated in CFU-E cells on trajectory A, as anticipated (Figure 6A-B; supplemental Tables 5-6). Heme-related genes, but not p53 pathway genes, were upregulated in trajectory-B cells of patients and controls (Figure 4B-C). Furthermore, of note, p53 activity is not increased in Flvcr1-deleted mice⁵⁹ in which excessive intracellular heme results from failed heme export. Together, this indicates that excess heme and not p53-dependent mechanisms triggers the demise of trajectory-B cells.

DBA is an inherited, germ line disease, and, as such, all DBA cells contain the mutation resulting in ribosomal protein haploinsufficiency. However, MDS-5q is an acquired clonal disease that originates in a hematopoietic stem/progenitor cell and results in deletions within chromosome 5q and haploinsufficiency of RPS14. One pathophysiological dilemma is that at presentation, patients with MDS-5q have anemia, even when many unmutated healthy hematopoietic stem/progenitor cells remain. At the time of sample collection, 74.5% of marrow cells in a patient with MDS-5q (patient M1) had a chromosome 5q^–, when assayed via interphase fluorescence in situ hybridization (iFISH). Although 25.5% of the patient’s marrow cells had an intact chromosome 5, the patient had a profound macrocytic anemia (hematocrit = 20.9%, hemoglobin = 7.0 g/dL, mean corpuscular volume = 117 fL; normal ranges, 38.3%-48.6%, 13.2-16.6 g/dL, and 80-100 fL, respectively). Using expression ratios of RPS14 to RPS18, we calculated that 77% of erythroid cells contained the 5q^– (Figure 7A-C; supplemental Figure 10). Studies on samples from the second patient with MDS-5q were similar, with 65% of marrow cells (via iFISH) and 63% of the erythroid cells (based on expression) having 5q^– (Figure 7A-C; supplemental Figure 10), confirming the accuracy of this methodology.

We, then, separately analyzed the transcriptomes of the cells on trajectory A with a 5q^– chromosome and cells with a 5q⁺ chromosome (Figure 7D-E). The transcriptomes of the 5q⁺ cells were not normal but showed changes also present in the 5q^– cells (and DBA samples). For example, both the 5q⁺ and 5q^– cells showed early upregulation of HIF1A targets, ROS response, cell cycle regulation, and spliceosome pathways and downregulation of erythroid differentiation and mammalian target of rapamycin signaling pathways (Figure 7D-E; supplemental Tables 7 and 8). The 5q⁺ cells uniquely upregulated oxidative phosphorylation–related pathways, which were downregulated in the 5q^– cells. Conversely, the 5q^– cells upregulated several p53-associated pathways whereas the 5q⁺ cells only upregulated 2 p53-related genes, CDKN1A and BAX. Their upregulation in the 5q⁺ early CFU-E cells, suggests that p53 activation might contribute to the failure of the 5q⁺ cells to expand and complete erythroid differentiation. When we looked specifically at the expression of heme-regulated genes,⁶¹ we found the 5q^– cells on trajectory A had significantly elevated expression of heme-responsive genes at the early CFU-E stage (Figure 7F). Trajectory-A 5q⁺ cells were able to compensate sufficiently that heme-responsive genes were not globally increased (Figure 7G); however, the upregulation of oxidative phosphorylation pathways^60,66 (Figure 7E) suggests that there may still be elevated heme levels.

When CFU-E transition to proerythroblasts, iron is plentiful, and heme synthesis is robust, yet globin is insufficient, and cells depend on FLVCR to export heme and avoid heme- and ROS-mediated damage.^3,13-15 When FLVCR is conditionally deleted in mice, the animals develop a severe macrocytic anemia characterized by excess heme in erythroblasts, excess cytoplasmic ROS, and cell death via apoptosis and ferroptosis.^14-17 Previously, we,³ and others,^63,67,68 hypothesized that the macrocytic anemia of DBA and MDS-5q, disorders associated with germ line or somatic ribosomal protein haploinsufficiencies, respectively, have a related pathogenesis. Here, using a combined data set to minimize the impact of individual variability, we explored the common pathophysiology of these disorders with single-cell analyses. The single-cell pseudotime analyses allowed us to dissect the decision-making of individual cells at designated points along their differentiation.

Uncultured cells (day 0) from patients with DBA and MDS-5q and cells after 3 and 6 days of culture distributed onto 2 trajectories, a healthy trajectory (trajectory A), resulting in effective erythroid differentiation, and an unhealthy trajectory (trajectory B), resulting in ineffective erythroid differentiation and cell death. Data analyses and succinylacetone treatment (Figures 5-6) confirmed that heme excess drives a cell toward trajectory B and its ultimate demise. Interestingly, some (24% ± 12%, range, 14%-44%) erythroid marrow cells from healthy controls also follow trajectory B (supplemental Table 4) and have gene expression patterns just like patient trajectory-B cells (Figure 4; supplemental Table 3), suggesting that heme is a broad regulator of effective vs ineffective erythropoiesis. The reason why some cells from healthy individuals follow trajectory B is uncertain, but this potentially assures a ready supply of early progenitors so that erythropoiesis can rapidly expand at times of physiological need.

By the C2 stage, the cells on trajectory B have significant upregulation of heme-regulated genes (Figure 5D) but fail to upregulate erythroid differentiation genes, indicating that a primary cause of death is excess heme. That reducing heme synthesis improves cell survival (Figure 5E) suggests that excessive heme is a key driver of apoptotic and ferroptotic cell death and that heme/globin balance at each stage is essential for effective erythroid differentiation. There is substantial time between when cells begin to follow trajectory B and when they die, providing a window for therapeutic intervention, for example with drugs that diminish heme synthesis or chelate iron.

Cells from patients with DBA or MDS-5q that follow the healthy trajectory A are not normal. Rather, they have successfully compensated to avoid fatal damage. Firstly, they have reduced erythroid differentiation genes (including heme biosynthesis genes), which is a direct response to excess heme.^59,60 Secondly, they express increased DNA damage response genes, likely to facilitate repair of ROS-mediated cell damage. And thirdly, they have upregulated HIF1A targets and hypoxia response genes, which may facilitate cell survival in light of the reduced heme synthesis. Many of these findings have been previously validated both in vivo and in vitro.^3,59,60 That culture stabilizes precarious DBA and MDS-5q cells allowing a greater number to survive and effectively differentiate (follow trajectory A) suggests that factors present in culture help cells compensate. That some patients with DBA spontaneously remit⁶⁹ implies that cell-extrinsic factors can also improve the in vivo differentiation of erythroid cells haploinsufficient for a ribosomal protein. It seems possible that drugs that mimic these compensations might ameliorate DBA or MDS-5q anemia.

We also found that the cells with 5q⁺ present in patients with MDS-5q have significant alterations in gene expression despite normal levels of RPS14 (Figure 7E; supplemental Table 8). The commonly upregulated pathways (ROS response and cell cycle regulation pathways) suggest that both 5q⁺ and 5q^– cells are responding to oxidative damage. This would indicate that the 5q^– cells damage the marrow microenvironment and indirectly impair the differentiation of neighboring genetically normal 5q⁺ cells, blocking their optimal expansion. The 5q⁺ cells partially compensate but remain affected and, thus, cannot fully expand and differentiate (Figure 7E,G).

Our studies also provide insight into the posttranscriptional mechanisms that regulate erythropoiesis, such as the alternative splicing of mRNA.^51-53 In addition to mature spliced transcripts, there was stage-specific and gene-specific intron retention (Figure 3B-C), leading to transcripts unlikely to translate into proteins. There was such a prevalence of intron retention that RNA velocity analysis predicted a reversed trajectory of differentiation after the CFU-E stage (Figure 3A). Besides providing a vehicle to order cells and align studies, the data imply that alternative splicing and intron retention may be especially important regulatory processes in red blood cell differentiation, as proposed by others.^45,46,51-53 

The authors thank the patients, especially patient D1 for their support of this project. The authors thank Chenkai Luo for early work developing and testing demultiplexing methods, and Kalyanashis Chakraborty for preliminary analysis of external data sets and integrating them into our analysis pipeline.

This research was funded by National Institutes of Health grants HL031823 (J.L.A.) and CA190122 (Q.T.), and the Institute for Systems Biology’s Innovator Award Program (C.G.L.). Human bone marrow samples were obtained at the University of Washington clinical research center, which is supported by National Institutes of Health National Center for Advancing Translational Sciences award number UL1TR000423.

Contribution: R.T.D., J.L.A., and Q.T. designed the studies; R.T.D., A.D.M., Z.Y., and C.G.L. performed the experiments and generated figures; C.G.L. processed the raw data and developed novel analytical methods; C.G.L., C.M., and X.Y. performed the bioinformatics analysis; R.T.D. and J.L.A. coordinated and directed the studies; R.T.D., C.G.L., and J.L.A. wrote the manuscript; and all authors edited the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

The current affiliation for Z.Y. is World Laureates Association (WLA), WLA Laboratories, Pudong New District, Shanghai, China.

The current affiliation for Q.T. is ProfoundBio, Seattle, WA.

Correspondence: Janis L. Abkowitz, Division of Hematology, Department of Medicine, University of Washington, Box 357710, 1959 NE Pacific St, Seattle, WA 98195; e-mail: janabk@uw.edu.