The functional landscape of alternative splicing in hematopoietic lineage commitment

Domain: neuroscience
Published: 2026-01-27
Editorial outcome: accept

Abstract
Alternative splicing (AS) is a ubiquitous post-transcriptional regulatory mechanism, that has greatly expanded the transcriptomic and proteomic diversity in vertebrates. While gene regulation of hematopoiesis has been extensively researched in vertebrates, the functions of species- and cell lineage-specific splice variants in vertebrates are largely unknown. Here, we curate transcriptomic data on fetal hematopoietic organ development in six vertebrates and hematopoietic cell differentiation in humans and mice. To identify functional exon-skipping events among thousands of cassette exons in protein-coding genes for a specific differentiation lineage and species, we develop a machine-learning model interrogating 19 features including dynamic expression, protein structure, and evolutionary conservation, and integrate them into a single prediction score, named Functional AS Score (FAScore). Using FAScore, we identify four previously-uncharacterized functional AS events in which deletion of the AS exon leads to defects in erythropoiesis and myelopoiesis. Furthermore, we demonstrate that deletion of exon 15 of TBC1D23 reduces erythropoiesis in mice and zebrafish through elevated binding capacity to RANBP2/RANGAP1 leading to increased SUMOylation level of HDAC1. Collectively, our study presents a valuable tool to identify functional exon skipping (ES) events during hematopoietic lineage commitment, and establishes a research paradigm that can be broadly applied to other biological processes.

Introduction
Alternative splicing (AS), a key post-transcriptional regulatory mechanism, increases the complexity of eukaryotic transcriptomes and proteomes. High-throughput transcriptomic analyses unveiled that 95% of multi-exonic genes undergo AS in human1. AS not only contributes to transcript diversity but also plays a critical role in regulating essential physiological processes, including organ development including brain2–4 and sexual organs5,6, as well as cell differentiation processes2–9 such as hematopoiesis10–13. In line with this, abnormal AS can influence disease states14–16 and consequently understanding the role of AS is essential for developing gene repair-based therapeutic strategies.

Hematopoiesis has been extensively studied as a model for understanding stem cell biology and developmental processes. In mammals, hematopoiesis begins during the embryonic stage, with hematopoietic stem cells (HSCs) first emerging in the aorta-gonad-mesonephros (AGM) region. As development proceeds, these cells migrate to the fetal liver (FL), and then to bone marrow (BM). Hematopoiesis is the process by which various blood cell lineages (collectively referred to as ‘lineage’ thereafter) are produced. We categorize the mature blood cells derived from the same progenitor cells into three lineages, including erythroid, myeloid and lymphoid lineage. The erythroid lineage, originating from megakaryocyte-erythrocyte progenitor (MEP), encompasses erythrocytes and megakaryocytes. The myeloid lineage, arising from granulocyte-monocyte progenitor (GMP), comprises morphologically phenotypically, and functionally distinct cell types, such as granulocytes, macrophages, and monocytes. The lymphoid lineage, derived from common lymphoid progenitor (CLP), includes T cells, B cells, and natural killer (NK) cells. AS plays significant roles in regulating the functional diversity of hematopoietic cells in human and mouse, including HSC formation7,9,17, erythropoiesis18,19, megakaryopoiesis19,20 and granulopoiesis21. Despite extensive efforts, only a small fraction of the AS introduced transcript diversity has been functionally validated as it remains experimentally challenging.

The advancement of machine learning (ML) algorithms provided new solutions to predict the functions of isoforms22–29. However, previous ML approaches were mainly based on static genomic information, such as conservation, sequence motif, protein domain, structure or AS level in certain cell types, tissues, organs or diseases17,30,31, and lacked dynamic expression information of the isoforms. Cross-species dynamics for genes32, lncRNAs33, and AS34 has proven to be a powerful framework to assess their evolutionary conservation in mammalian organ development. Thus, incorporating data on both dynamic and static features may offer a better strategy to screen for functional AS. Here, we developed a machine-learning model called functional AS score (FAScore), interrogating 19 features, including dynamic expression pattern, evolutionary conservation, and protein structural information. Since ES events represents the most prevalent type of AS and has been reported to contribute more functional events compared to other types, we focused our model on this AS type. FAScore method used (a) PU-learning coupled with rebalancing methods to address the challenges of limited positive sets, absence of negative sets, and data imbalance, and (b) Random Forest (RF) to screen for functional AS involved in organ development and hematopoietic differentiation. Each AS event has a single FAScore of functional relevance in a given differentiation lineage and species. As an example of the application of FAScore, we identified four AS events that were important for specific hematopoietic lineage commitment.

Results
Dynamically regulated alternative splicing atlases of vertebrate hematopoiesis
To study the evolution of alternative splicing of hematopoietic lineage commitment, we integrated 271 bulk RNA-seq samples spanning hematopoietic stem and progenitor cells (HSPCs) and fetal hematopoietic organ (FHO) development. The HSPC dataset includes 17 cell types belonging to three lineages in mouse bone marrow (BM)35 and human umbilical cord blood17 (Supplementary Data 1). We also analyzed datasets of the fetal liver, the major fetal hematopoietic organ36 across five mammals (human, rhesus, rabbit, rat and mouse)32. In addition, we included the caudal hematopoietic tissue (CHT) which is the embryonic hematopoietic organ in zebrafish from 36 h post fertilization (hpf) to 4 days post fertilization (dpf)37 (Supplementary Data 1, and Fig. 1a).

We identified and quantified alternatively spliced exons and AS events as described in MeDAS (Metazoan Developmental Alternative Splicing database)38. The percent spliced-in index (PSI) for each exon or ‘exon part’ (exon-centric PSI, ePSI) was calculated. AS types were subsequently inferred using SUPPA239 by comparing exon regions with reference annotations (see “Methods”). Additionally, we identified 1,640,242 splicing junctions (SJs), of which ~11.9% were unannotated canonical SJs (Supplementary Fig. 1a). Exon skipping (35.7%) was the most frequent AS types (Supplementary Fig. 1b). On average, 17,110 cassette exons (from 2953 to 32,400), originating from 6390 protein-coding genes (from 1824 to 9015), were contained across species and lineages (Supplementary Data 2). Given the large numbers of cassette exons and more functional ES events were documented than other AS types, we focused on ES in the subsequent analyses.

Genes or AS events that are dynamically regulated during developmental or cell differentiation trajectories are more likely to play critical roles in these biological processes, as their temporal regulation often reflects direct functional involvement in key cellular transitions or fate decisions34. To quantify the dynamic changes in gene expression and alternative splicing, we computed a series of feature scores based on the expression levels or PSI values across developmental timepoints or differentiated cell types (Fig. 1b). Specifically, we calculated the expression breath (Gene.Range and AS.Range) to assess the magnitude of change, the correlation coefficient (Gene.cor and AS.cor) and its statistical significance (Gene.cor.p and AS.cor.p) to evaluate the association with temporal or differentiation trajectories, the linear regression slope (Gene.slope and AS.slope) and its statistical significance (Gene.slope.p and AS.slope.p) to determine the rate and the statistical significance of the progressive changes, as well as tau (Gene.tau and AS.tau) to measure expression specificity (Fig. 1b, see “Methods”). As a representative example, the MEIS1, a transcription factor whose inactivated impaired hematopoietic stem cell niche development, megakaryocyte and platelet deficiencies in mice40–42. We quantified six features of the dynamic gene expression of MEIS1 in human erythroid lineage. It showed a gradually decreased expression in hematopoiesis with an upregulation in megakaryocyte (Fig. 1c), consistent with a previous report43. Meanwhile, MEIS1 contains two ES events, the exon (E) 1 and 12. The MEIS1-E12 exhibits high level of changes (Fig. 1d), consistent with its function in inducing the fate decision of human hematopoietic progenitors toward megakaryocyte-erythroid progenitors (MEPs)44. In contrast, the MEIS1-E1 showed a stable expression across cell types with low values of dynamic features (Fig. 1e).

To facilitate intuitive assessment of the dynamics of genes and AS events, we integrated the dynamic features described above into a unified dynamism score (DyScore). Based on its bimodal distribution, we applied Gaussian Mixture Modeling (GMM) to classified the genes or AS events into “dynamic” and “non-dynamic” (Supplementary Fig. 1c, e, see “Methods”). Overall, 27% of AS events and 56% of genes were classified as dynamic (Supplementary Fig. 1d, f). Notably, in both human and mouse, dynamically regulated AS events were more likely to be cell lineage-specific, with 63% dynamic in only one hematopoietic lineage compared to 29% of dynamically expressed genes (Supplementary Fig. 1g). Moreover, dynamically regulated AS exhibited significantly higher evolutionary conservation scores in the regions flanking the splice sites across all lineages and fetal hematopoietic organ (paired t-test, p<1.5×10-4, Supplementary Fig. 1h). These results suggest that dynamically regulated AS events may play specialized and conserved regulatory roles in hematopoietic lineage, highlighting the importance of splicing dynamics in the spatio-temporal control of transcriptome complexity during hematopoiesis.

Predictive machine learning model for functional alternative splicing
To screen for functional AS in hematopoiesis, we extracted features other than transcriptomic dynamics including sequence conservation, structural homologs and integrity, domain integrity, transmembrane helices, signal peptide and target peptide. Using these features, we developed a machine-learning model named FAScore (Functional AS Score) to predict functional AS events in a given differentiation lineage and species (Fig. 1a). For each AS event, the model assigns a species- and lineage-specific FAScore, designated as species_lineage_FAScore (e.g., Human_erythroid_FAScore, Mouse_erythroid_FAScore). Due to the lack of a verified negative set and the limited number of known functional AS events with experimental evidence (37 unique exons, with redundancies due to cross-species conservation and multi-lineage functionality, Supplementary Data 3, see Methods), we employed a Positive-Unlabeled (PU) Learning framework to identify reliable positive and negative sets for model training (see Methods). It defined a positive set (4803) and a negative set (5245) to form the training dataset (Supplementary Fig. 2a). A comparative analysis of feature scores among the gold-standard positives, the reliable positives, and the reliable negatives revealed minimal differences between the gold-standard and reliable positive sets, but significant differences between the reliable positives and negatives, supporting the effectiveness of PU-learning (Supplementary Fig. 2b). Compared to Naive Bayes (AUC: 0.94), Logistic Regression (0.96), Support Vector Machine (SVM, 0.98), AdaBoost (0.98) and CART (0.90), the Random Forest (RF) model demonstrated the best performance (0.99). Therefore, we developed a RF model on this constructed training set and selected the best-performing model by 10-fold cross-validation (AUC-ROC=0.993, AUC-PR=0.992, F1 score=0.952, Rank index=0.243, FNR=0.050) (Fig. 2a, and Supplementary Fig. 2c). We assessed the contribution of each feature to the prediction by calculating the relative importance score using Gini coefficient45. The features contributing the most to our model were dynamic splicing (AS.slope.p, 21% and AS.cor.p, 12%), sequence conservation cross-species (14%) and protein domain (Domain.integrity, 10%) (Supplementary Fig. 2d). Our model integrates multiple features into a unified predictive score, offering a comprehensive selection criterion to screen for functional AS which may not be identified focusing on features of conservation or protein domain (Supplementary Fig. 3a). For example, using conservation as the sole criterion would exclude low-conservation yet functional AS events found in the known positives set, such as Kdm1a-E746 and EHMT2-E1047 (Supplementary Fig. 3b). Moreover, the conservation and FAScore are weakly related as there are events with high conservation but having low FAScore and events with high FAScore with low conservation (Supplementary Fig. 3c, d).

Next, we measured the FAScores for annotated and unannotated AS events from fetal hematopoietic organ development and hematopoietic cell differentiation per species (Supplementary Data 2), that were not used in the model training deriving a continuous score for AS, ranging from 0 to 1, with 1 as most likely to be functional. By fitting a three-component GMM, we defined the following classes by probability density intersection method: ‘functional AS’ (FuncAS), ‘non-functional AS’ (NonfuncAS), and ‘uncertain AS’ (Uncertain) (Supplementary Fig. 4a). In total, we predicted 30,615 functional AS (16.3%), 53,454 non-functional (55.3%), and 104,133 uncertain (28.4%) (Supplementary Fig. 4b, and Supplementary Data 2). Furthermore, we compared the 19 features among these three groups and observed significant differences in features between the predicted functional and non-functional sets, consistent with the observation in the model training (Supplementary Fig. 4c).

We propose that AS events with highly conserved splicing sites and harboring functional domains or modification sites are more likely to play a role in complex biological processes such as hematopoiesis. To evaluate the prediction outcomes, we compared the predicted functional set (FuncAS) and non-functional set (NonfuncAS) from four distinct aspects, including protein domain, conservation, modification sites and functional pathways of harboring genes. First, we observed that compared to the non-functional, the functional AS were more likely to be annotated with functional domains in the UniProt database48 such as transmembrane domain, active sites, protein-protein interaction, intrinsic disorder region (prop.test, p<0.05) (Fig. 2b, and Supplementary Fig. 5a). Second, the sequences around the splicing junctions were more conserved in functional AS than in the non-functional AS events assessed by PhastCons score34 (Fig. 2c, and Supplementary Fig. 5b). Third, functional AS were more likely to harbor post translational modification (PTM) sites annotated in PhosphoSitePlus49, such as methylation, phosphorylation and ubiquitination (prop.test, p<0.05) (Fig. 2d, and Supplementary Fig. 5c). Finally, we found the genes that harbour at least one functional AS event (FuncAS-harbouring genes) tended to be enriched in splicing factors (SF), transcription factors (TF) and RNA binding proteins (RBP) (Supplementary Fig. 5d). GO analysis of these genes are associated with pathways related to the lineage, as ‘erythroid differentiation’ was prominently enriched in the erythroid lineage; ‘regulation of lymphocyte activation’ was enriched in the lymphoid lineage; and ‘myeloid cell differentiation’ along with ‘neutrophil activation’ were enriched in the myeloid lineage (Fig. 2e). In particular, for fetal hematopoietic organ development, the FuncAS-harbouring genes were enriched in pathways related to both liver development and hematopoiesis. Geneset enrichment analysis (GSEA) also showed that FuncAS-harbouring genes were enriched with transcription factors (TFs) of the myeloid, such as KLF4, IRF8, STAT1 and STAT3 (Supplementary Data 4, and Supplementary Fig. 5e). These results suggest the accuracy of our prediction model.

FAScore predicts cell lineage-specific AS events
To identify conserved AS, we next mapped the genomic locations of 42,452 AS events hosted in 11,075 genes identified in the human dataset and found 16,627 AS events had at least one corresponding AS region in other vertebrate species. These AS events were then classified as primate-specific (age class 1), mammal-specific (class 4) or ancient (present in all vertebrates, class 5, dating back 429 million years ago (MYA) (Fig. 3a).

The creation of novel AS exons occurred fairly often during vertebrate evolution. For example, the 1,833 primates-specific AS events were only present in humans and macaques but not in other species (Fig. 3a). Using FAScores, we observed a marked increase in the proportion of predicted functional AS events in class 5 compared with the other classes (Fig. 3b), but not in non-functional AS exons. In total, 6,748 AS events from 3,073 genes were identified as conserved in vertebrates (class 5). It is reported that older genes have more critical functions50, and our results suggest that ancient AS events may also play essential roles (Fig. 3b).

Hematopoiesis in the vertebrates takes place at various times and anatomical locations, including fetal liver/CHT (defined as FHO) and BM. Analyzing the spatiotemporal expression patterns of the conserved AS, we identified 24 clusters from the 6,748 vertebrate-conserved AS events (age class 5) using Mfuzz51 (Supplementary Fig. 6a, and Supplementary Data 5). About 27.2% of AS events were functional exclusively in FHO (C2=1,837), while 8.8% (C6, C12 and C17, 597) were functional in all three lineages postnatally in BM (Supplementary Fig. 6a). Meanwhile, 40.1% AS events were cell lineage-specific in hematopoiesis, including 19.4% myeloid-specific (C4, C5, C15 and C22), 9.6% erythroid-specific (C8, C11 and C21) and 11.1% lymphoid-specific (C9, C10 and C14) (Fig. 3c, Supplementary Fig. 6a). These findings indicated that key AS events are likely dynamically regulated during lineage commitment and in different anatomical locations.

Next, we selected candidates with top FAScores from different cell lineages of mouse or human for experimental validation (Fig. 3d). For the erythroid lineage, TBC1D23-E15 (Human_erythroid_FAScore: 0.998) and EPB41L1-E14 (Mouse_erythroid_FAScore: 0.995) were chosen. For the myeloid lineage, KLF6-E3 (Human_myeloid_FAScore: 0.997) and SSBP3-E6 (Mouse_myeloid_FAScore: 0.915) were chosen. And for the lymphoid lineage, FYN-E7 (Mouse_lymphoid_FAScore: 0.987) was selected (Supplementary Fig. 6b). Given that FYN-E7 has been reported to play an important role in T cell signal transduction in mice52, the other four events were experimentally analyzed. SSBP3-E6 and KLF6-E3 are more frequently included in the myeloid lineage (average PSI: 0.27 and 0.92, respectively) compared to the erythroid (average PSI: 0.08 and 0.77) and lymphoid (average PSI: 0.10 and 0.75) lineages (Fig. 3e). In contrast, TBC1D23-E15 is more likely to be included in the erythroid lineage (average PSI: 0.78), while tends to be skipped in the lymphoid (average PSI: 0.32) and myeloid (average PSI: 0.39) lineages (Fig. 3e). We first validated that all of the candidate exons were all transcribed in the colony-forming unit-granulocyte-macrophage (GM), the colony-forming unit-granulocyte-erythrocyte-monocyte-megakaryocyte (GEMM) and the burst-forming unit-erythrocyte (BFU-E) colonies derived from colony-forming unit (CFU) assays (Supplementary Fig. 6c). To validate that they are functional in given lineages for which they have high FAScores, we performed CFU assays on the mouse HSPCs in which the candidate exons were deleted using GFP-expressing retroviral vector having two sgRNAs targeting the flanking intronic regions of the exons (Fig. 3f). The transduced cells with GFP expression were sorted and plated for CFU assays. The CFU assays revealed that the deletion of SSBP3-E6 or KLF6-E3 in mouse c-kit+ cells resulted in a reduced formation of myeloid colonies by 58% and 38% respectively, while the erythroid colonies were of wildtype level (Fig. 3g). Similarly, the CFU assays showed that the deletion of the Tbc1d23-E15 inhibited erythroid colony formation by 35% without a significant impact on myeloid populations (Fig. 3g). EPB41L1-E14 which exhibits a higher inclusion rate in mature myeloid cells (Fig. 3e) and the deletion of Epb41l1-E14 markedly enhanced erythropoiesis with an increase of 38% of BFU-E colonies while no significant difference on the myeloid colonies was observed (Fig. 3g). In summary, FAScore accurately identified AS events play important roles in lineage fate commitment.

Loss of TBC1D23-E15 inhibits vertebrate erythropoiesis
As shown above, FAScore aids in identification of ES events that may play important roles in hematopoiesis, we next experimentally analyzed the role of TBC1D23-E15 with high Human_erythroid_FAScore (0.998) and high PSI (0.78) in erythroid lineage of human (Fig. 3d, e). TBC1D23, a member of the TBC family, is known to regulate participate in the regulation of the traffic between endosomes and Golgi apparatus53,54 and is reported to influence the pathway of the innate immune response55 (Supplementary Fig. 7a). We characterized the transcription level of Tbc1d23-E15 in the sorted HSPCs and mature cell types in mouse bone marrow using RT-PCR, the result show the PSI correlated well with the PSI derived from bulk RNA-seq analysis of mouse hematopoietic cell populations, confirming that E15 exhibits the most pronounced differential splicing in erythroid cells (Supplementary Figs. 7b, c). By comparing the adjacent regions of TBC1D23-E15 in 13 species ranging from C. elegans to humans, we found that the canonical splice site dinucleotide motifs (AG at the 3’ splice site and GT/GC at the 5’ splice site) were present in zebrafish but not in amphioxus, suggesting that this exon arose during the evolutionary transition from invertebrates to vertebrates (Supplementary Fig. 7d). Using the human erythroleukemia K562 cell line and TBC1D23 antibody, we conducted a mass spectrometry (IP-MS) experiment and confirmed that this 45-bp exon indeed encodes a peptide of 15 amino acids (Supplementary Fig. 7e). This peptide is proximal to the binding region with FAM91A153,56.

To investigate the role of TBC1D23-E15 in vertebrate hematopoiesis, we generated a transgenic mouse line in which Tbc1d23-E15 was deleted (Fig. 4a). In contrast to wild type (WT) mice expressing both long and short isoforms of Tbc1d23, the Tbc1d23ΔE15/ΔE15 (E15-KO) mice express only the short isoform (Fig. 4b, and Supplementary Fig. 7f) with the overall Tbc1d23 protein level was comparable to WT mice (Fig. 4c, d). E15-KO mice exhibited significant reductions in red blood cell (RBC) numbers, hematocrit (HCT) and hemoglobin (HGB) levels compared to the WT littermate controls while white blood cell counts remained unchanged (Fig. 4e, f). In addition, to assess the conservation of this phenotype, we targeted the long isoform of Tbc1d23 in zebrafish (Tbc1d23-L, including E15) by morpholino (MO) treatment. The Tbc1d23-L-MO (E15-MO) treatment specifically reduced the mRNA level of the long isoform while increasing the expression of the short isoform, with no significant effect on total Tbc1d23 protein level comparable to the fish treated with the control MO (Ctrl-MO) (Fig. 4g–i, Supplementary Fig. 7g). Using transgenic fish expressing DsRed under the control of the Gata1 promoter, we found a significant (45%) reduction in RBC levels in 48 hpf fish under E15-MO treatment, compared to control (Fig. 4j, k). Consistent with this, O-dianisidine staining revealed an early 50% reduction of HGB levels in E15-MO-treated fish (Fig. 4l, m).

We next determined whether Tbc1d23-E15 deletion leads to a defect in erythroid differentiation in BM which results in anemia in mice. Firstly, we found that the number of bone marrow cells (BMCs) in the E15-KO mice was 52% lower than in WT mice (Fig. 5a). Moreover, by characterizing the stages of erythroid cells (referred to here as RII, RIII, RIV and RV) by flow cytometry utilizing TER119 and CD44, we observed the frequency of erythroid cells (RII, RIII, RIV and RV) in BM was also halved compared to WT mice (Fig. 5b). In contrast, the proportion of HSCs and progenitors in the BM of E15-KO mice remained unchanged (Supplementary Fig. 7h). Colony formation assay showed a 23% reduction in BFU-E colony number and a 40% reduction in colony-forming unit-erythroid (CFU-E) clones for E15-KO mice (Fig. 5c, d).

To further investigate how Tbc1d23-E15 affects erythropoiesis, we performed a single-cell RNA sequencing (scRNA-seq) analysis on BMCs of E15-KO and WT mice (Fig. 5e, and Supplementary Fig. 8a). In total, we identified 802 differentially expressed genes (DEGs) in all cell types, with the greatest number of DEGs observed in the pro-erythroblasts (199) (Supplementary Fig. 8b, c). To identify the lineage with the most transcriptional changes between E15-KO and WT, we used Augur57 to calculate AUC values. Our analysis revealed that erythroid cells displayed a higher AUC value than other lineages, indicating a strong response to E15-KO (Supplementary Fig. 8d). Utilizing CytoTRACE58, we assessed the differentiation potential scores of all cell types and found that pro-erythroblasts received the highest score (0.46), compared to 0.10 of MEP (Fig. 5f), in line with our observations that deletion of Tbc1d23-E15 affects the maturation of erythroid cells but not progenitors.

TBC1D23-E15 modulates erythropoiesis by regulating the SUMOylation of HDAC1
To investigate how TBC1D23-E15 regulates erythropoiesis, we analyzed the interactomes of long and short isoforms of TBC1D23 by an integrated multi-omics approach. We generated three in both K562 and HEK293T, one lacking TBC1D23 expression (KO) and the other two expressing only TBC1D23-L or TBC1D23-S (Supplementary Fig. 9a, and Fig. 6a, b), and performed IP-MS in the K562 cell lines (Fig. 6c, Supplementary Data 6). RanBP2 was identified as one of the candidate proteins showing differential binding toward TBC1D23-S or TBC1D23-L (Supplementary Fig. 9b, c). RanBP2 is a known SumoE3 ligase that forms a complex with SUMO1-modified RanGAP159. The immunofluorescence analysis showed that TBC1D23 co-localized with RanBP2 (Fig. 6d). Using the TBC1D23-KO, TBC1D23-L and TBC1D23-S K562 cell lines overexpressing RanBP2 or RanGAP1 E3 SUMO ligase domain, we confirmed TBC1D23-S consistently precipitated more RanBP2 and RanGAP1 than TBC1D23-L by 5 and 4 times, respectively (Fig. 6e, f), suggesting that E15 exclusion enhances engagement with the complex (Fig. 6e, f).

To investigate whether the SUMOylation level was differentially regulated by the long and short isoforms of TBC1D23, we performed quantitative SUMO proteomics using K562 cells expressing only TBC1D23-S or -L (Fig. 7a, and Supplementary Data 7) and TBC1D23-KO cell lines were also included as a negative control (Supplementary Data 7). The identified substrate proteins were enriched in various biological processes, with ‘regulation of hemopoiesis’ ranked as one of the top 10 pathways (p=5.2e-04) (Fig. 7b). Eleven proteins within this pathway showed differential SUMOylation in TBC1D23-S and -L expressing cells, and HDAC1, TFE3 and HMGB1/HMGB2 had higher SUMOylation levels in cells expressing TBC1D23-S compared to those expressing TBC1D23-L (Fig. 7c). To elucidate the potential key regulators, we further performed consensus co-expression network analysis on the scRNA-seq dataset using hdWGCNA60 (Supplementary Fig. 10a). The most significantly changed network out of the seven identified co-expression modules (M) were M1 and M5 (Supplementary Fig. 10b, c). M1 was associated with erythrocyte differentiation and protein modification, while M5 was associated with immune response (Supplementary Fig. 10d). Notably, HDAC1 (histone deacetylase 1) and PRMT1 were present both in M1 and the differential SUMOylation dataset related to the regulation of hemopoiesis (Supplementary Fig. 10e).

Consistent with the tighter binding between TBC1D23-S and the RanBP2/RanGAP1 complex, the SUMO1-mediated SUMOylation level of HDAC1 was 5 times higher in cells expressing TBC1D23-S than those of TBC1D23-L. Importantly, TBC1D23-S and -L displayed similar binding toward HDAC1 (Fig. 7d, e, Supplementary Fig. 10f). To confirm that HDAC1 is a substrate of RanBP2, we expressed HDAC1 with RanBP2 wt or a catalytic-dead mutant, RanBP2mut (F2735A/F2736A/C2737A), and found that HDAC1 was drastically SUMOylated only in RanBP2wt-expressing cells (Supplementary Fig. 10g, and Fig. 7f). Lys476 of HDAC1 was identified in our assay with higher SUMOylated levels in cells expressing TBC1D23-S. Lys476 was reported as a major SUMOylation site of HDAC1, together with Lys44461. To determine whether the two sites were SUMOylated by the RanBP2/RanGAP1 complex, we generated cell lines expressing FLAG-tagged HDAC1 wt (K444/K476) or HDAC1 mut (R444/R476) and found that RanBP2-dependent SUMO1-mediated SUMOylation can only be detected in HDAC1 wt expressing cells (Fig. 7f). Furthermore, HDAC1 SUMOylation was elevated by 1.5-fold in BMCs from E15-KO mice compared to those of WT controls (Fig. 7g, h). Importantly, the reduced CFU-E colony numbers of E15-KO BMCs could only be rescued by the expression of HDAC1 mut, but not by that of HDAC1 wt (Fig. 7i). Thus, TBC1D23-E15 regulates erythropoiesis likely by modulating the SUMOylation of HDAC1.

To reconstruct the gene regulatory networks (GRN) affected by Tbc1d23-E15, we analyzed the activity of transcription factors (TFs) in erythroid scRNA-seq data of E15-KO and WT BMCs by SCENIC62 (Supplementary Data 8). This analysis yielded 204 regulons out of 1,721 initial co-expression modules, with significantly enriched motifs for the corresponding TFs in erythropoiesis. Twelve TFs were significantly overlapped between these TF regulons and Hdac1-related TFs of erythropoiesis63 (p=0.003, Supplementary Fig. 10h, and Supplementary Data 9). For example, the downregulation of Myc is essential for terminal erythroid maturation64, whose activity increases in E15-KO and is associated with 2,774 predicted target genes. In contrast, Klf1, a master erythroid gene regulator65, shows decreased activity in E15-KO, with 1,198 predicted target genes (Supplementary Fig. 10i, and Supplementary Data 9). Furthermore, we found five correlation networks amongst different TFs by hierarchical clustering, in which TFs targeted by Hdac1 mainly belonged to network 1 (N1) (Spi1, Gata2, Myc, Stat5a, Rcor1) and network 2 (N2) (Tal1, Klf11, Klf1, Gata1, Sp1) (Supplementary Fig. 10j). Additionally, we obtained a co-regulation network conducted by the 10 TFs and their target genes to affect erythroid differentiation (Supplementary Fig. 10k). These results indicate that skipping E15 of Tbc1d23 leads to SUMO-modification of Hdac1, which may modulate the deacetylation and activation of the target TFs, impairing erythroid cell differentiation. This suggests the significant role of this AS event in regulating transcriptional networks during erythroid lineage commitment.

Methods
FAScore algorithm overview
The framework of the FAScore consists of several steps: (1) calculation of the dynamic features of genes or alternative splicing events; (2) extraction of structural features of isoforms from the APPRIS database; (3) Training of predictive model; and (4) calculation of the probability to generate the functional AS scores (FAScore). We implemented and tested FAScore in R (v4.0) and designed it for an R package. FAScore code is open source and available on GitHub (https://github.com/LuChenLab/FAScore.git), along with detailed descriptions of functions and tutorials.

Step 1: Calculation of the dynamic features
Data collection. To train the functional AS predictor and study the evolution of lineage commitment and developmental AS, we leveraged some bulk RNA-seq datasets comprising 1,985 libraries spanning the hematopoietic cell differentiation (87) and organ development (1,898). The hematopoietic differentiation datasets include 17 cell types belonging to three lineages extracted from the bone marrow (mouse)35 and umbilical cord blood (human)17 (Supplementary Data 1). And the organ development datasets were downloaded from Henrik Kaessmann et al.32 contain seven organs (forebrain/cerebrum, hindbrain/cerebellum, heart, kidney, liver, ovary and testis) from early organogenesis (mid-organogenesis for the heart) to adulthood across six mammals and chicken. In addition, we also exacted the dataset of the liver in the embryonic period across five mammals (human, rhesus, mouse, rat and rabbit) from the above datasets, and added the dataset of the caudal hematopoietic tissue (CHT) which was the embryonic hematopoietic organ in zebrafish from 36 hpf to 4 dpf stage37. Fetal liver and CHT are collectively referred to as the fetal hematopoietic organ (FHO) in the remainder of the article (Supplementary Data 1).

RNA-seq data preprocessing. All RNA sequencing libraries were processed using standard pipelines. Briefly, the quality of the raw sequence data was checked using FastQC software (v0.11.8). Low-quality reads were trimmed using Trimmomatic v0.3870 with parameters of ‘ILLUMINACLIP:adapters.fa:2:30:10 LEADING:3 TRAILING:3 and SLIDINGWINDOW:4:15 MINLEN:25’. Trimmed reads were aligned to the respective reference genomes (Supplementary Data 10) using STAR v2.6.1a71 with default parameters, supplemented by ‘-quantMode’ of TranscriptomeSAM, with the resulting BAM files, input into RSEM v1.3.172 to quantify the gene expression and isoform abundance as raw counts and transcripts per million (TPM) (Supplementary Fig. 11a–c). The strand specificity of each library was passed to the ‘-strandedness’ option of RSEM. The AS events were identified and quantified based on the exonic-centric method using the reported pipeline from MeDAS38 which identifies the exonic parts as AS events (Codes are available at https://github.com/LuChenLab/MeDAS). Constitutive exons exhibit a PSI of 1, while alternatively spliced exons range from 0 < PSI < 1. Based on their PSI values across all samples, exonic parts were categorized into: ‘not expressed’ (no PSI data available), ‘low inclusion’ (PSI < 0.05), ‘alternatively spliced’ (0.05 ≤ PSI ≤ 0.95), or ‘constitutive’ (PSI > 0.95). For downstream analysis, AS events were selected based on two criteria: (1) the exonic part must be classified as ‘alternatively spliced’ in at least two samples, and (2) the exonic part length must be ≥ 3 bp. Besides, AS types were inferred using SUPPA2 with ‘local AS events’ mode39, including ES, retained intron (RI), alternative 5’ splice site (A5), alternative 3’ splice site (A3), mutually exclusive exons (MX), alternative first exon (AF, resulting in mRNA isoforms with distinct 5’ UTRs), alternative last exon (AL). NA indicated novel AS events without annotated. For downstream analysis, we selected AS events whose annotated AS type included “ES” or those that were unannotated (novel).

Identification of dynamically expressed genes and AS events
Feature extraction. We extracted 6 features from the expression matrices, including the expression breadth (the range value, R), the expression relationship with differentiation or development (the Spearman’s rank coefficient, ρ; the Spearman’s test p-value, pρ), the expression specificity (the specificity index, τ) and the expression gradient (the linear regression coefficient, β; the p-value of the regression variable, pβ).
(1) Expression breadth: range value (R)
The value range (R) represents the difference between the maximum and minimum expression values, reflecting the extent of variation. For genes, expression breadth (Rgene) was calculated as follows: a pseudo-count of 1 was added to the TPM values, which were then log2-transformed; Rgene was then computed as the difference (max−min) among the transformed values for each gene (Supplementary Fig. 11d). when the expression breath is 0, it means no difference. Meanwhile, when expression breath is 1 or greater than 1, it means a minimum 2-fold expression difference (log2(FC)=1) within a gene. Thus, we applied a transformation step wherein any Rgene value ≥ 1 was set to 1 (Supplementary Fig. 11e). For AS events quantified using PSI, the range (RAS) was directly computed as the difference in PSI values for each AS event. Since PSI is naturally bounded between 0 and 1, RAS values also inherently fall within the 0–1 range (Supplementary Fig. 11f).
(2) Expression relationship with differentiation or development: Spearman’s rank coefficient (ρ) and Spearman’s test p-value (pρ).
The correlation of expression profile of development stages or cell types was calculated as genes or AS during organ development or hematopoietic stem cell differentiation. The developmental stages of organs or cell types along the lineage commitment were defined as ‘time-points’ starting from 1. To measure the relationship between the development or differentiation time-points and gene expression levels or AS PSI levels, we calculated two features, Spearman’s rank coefficient (ρ) and Spearman’s test p-value (Pρ, two-sided t-test)73. The rank coefficient ranges from −1 to 1, with the positive values indicating a direct correlation and the negative values indicating an inverse correlation. The p-values were transformed into binary form. The values below 0.05 were labeled as significant (1) and the rest as not significant (0).
(3) Expression gradient: regression coefficient (β) and p-value of the regression variable (Pβ)
To objectively characterize the linearity of relationships, we employed generalized additive modeling (GAM) and effective degrees of freedom (edf), which quantifies the complexity or “wiggliness” of a smooth term, reflecting the extent to which the model deviates from a straight line. The observation that edf values clustered near 174 indicates that gene expression and splicing levels exhibit predominantly linear tendencies (51.8%) (Supplementary Fig. 12a-c). Accordingly, we primarily conducted linear regression analyses, where the regression coefficient (β) represents the magnitude of influence exerted by the independent variable on the dependent variable, and the p-value of the regression variable (Pβ) reflects the statistical significance of this association.
To mitigate potential skewness in gene expression distributions while preserving biological interpretability, we applied a log2(TPM+1) transformation to gene expression values, which were then subjected to linear regression against normalized differentiation cell types or developmental time points. Although this transformation reduced bias from highly expressed genes, the resulting slopes initially ranged from −2 to 2, with most values concentrated near 0. These initial slopes showed only weak correlation with slopes derived from raw TPM values (Pearson r = 0.17, Supplementary Fig. 12d). We therefore applied capping to extreme values at the range [−1, 1], which affected fewer than 1% of genes. This adjustment substantially improved the correlation to r = 0.84 (Supplementary Fig. 12e), suggesting enhanced biological relevance. The optimized slope values also support direct additive integration across different features, providing a mathematically consistent basis for downstream dynamic score calculation. For AS events, PSI values inherently yield slopes within [−1, 1] without additional processing (Supplementary Fig. 12f), ensuring analytical consistency with processed gene expression metrics. Furthermore, p-values were binarized: values below 0.05 were assigned a significance label of 1, and all others were labeled 0.
(4) Expression specificity: specificity index (τ)
The specificity index (τ) was calculated to measure the expression specificity of gene/AS in organ development or lineage commitment75.
τ = (Σ (1−xi)) / (N−1) (1)
where N is the number of developmental stages or lineage commitment cell types, and xi is the gene expression or AS PSI normalized against the highest expression or PSI value in each gene or AS. τ interpolate the entire range between 0 for general and 1 for strict expression specific.

Definition of dynamically regulated genes and AS
The dynamic expression score (DyScore) was calculated as the arithmetic average of the six features described above. The DyScore is defined as:
DyScore = (R + |ρ| + Pρ + |β| + Pβ + τ) / n (2)
where R, Pρ, Pβ and τ represent the scaled range values, the Spearman’s test p-value, the p-value of the regression variable, and the specificity index, respectively. The |ρ| indicates the absolute value of the Spearman’s rank coefficient, and |β| is the absolute value of the regression coefficient. And n is the number of features. The DyScore value is between −1 to 1. The symbol indicates the direction of change, consistent with correlation.
Based on bimodal distribution of absolute dynamic score, we applied Gaussian Mixture Modeling (GMM) to classify the genes or AS events into “dynamic” and “non-dynamic”. A GMM with k=2 components effectively captured the bimodal distribution. The optimal threshold separating the two components was determined by probability density intersection method, resulting in a cut-off of 0.467 for genes and 0.324 for AS events. We then calculated the arithmetic mean of the gene-level and AS-level cutoffs, and it was 0.40. Finally, genes or AS events were divided into two classes, including dynamically expressed genes/AS events (≥0.4 or ≤ -0.4), and non-dynamically expressed genes/AS events (-0.4~0.4).

Step 2: Extraction of structural features
Structural characteristics of isoforms were derived from APPRIS database76 (https://appris.bioinfo.cnio.es), including 6 features:
(1) the structural homologs and integrity score. The presence of structural homologs in the Protein Data Bank (PDB) was evaluated and the integrity of the 3D structure was tested using Matador3D77.
(2) functionally important residues score. Conserved and functionally important amino acid residues were predicted by firestar78.
(3) domain integrity score. Domain integrity score was counted using the SPADEPfamscan program79.
(4) trans-membrane helices score. Transmembrane helices were predicted by three separate trans-membrane predictors of THUMP, including MEMSAT380, Phobius81 and PRODIV-TMHMM82 methods.
(5) signal peptide score. Signal peptides were predicted using SignalP.
(6) target peptide score. The presence of N-terminal presequences, including signal peptide (SP), mitochondrial transit peptide (mTP), chloroplast transit peptide (cTP) or thylakoid luminal transit peptide (luTP), were predicted using TargetP83.
The above features of annotations for five Ensembl species were downloaded, including human (GRCh38), rhesus (Mmul10), mouse (GRCm38), rat (Rnor6.0) and zebrafish (GRCz11).

Step 3: Training of predictive model
Benchmark dataset. AS events with actual evidence of functionality and phenotype were searched in PubMed84 using the terms of ‘alternative splicing’ with ‘brain’ OR ‘heart’ OR ‘liver’ OR ‘kidney’ OR ‘ovary’ OR ‘testis’ OR ‘hematopoiesis’ OR ‘myeloid’ OR ‘lymphoid’ OR ‘erythroid’. A total of 37 uniquely exon-skipping events were experimentally verified as functional in organ development or hematopoietic differentiation, among 2 AS events observed in two species and 15 observed in multiple lineages or organs (Supplementary Data 3), adding up to 58 positive instances. These were used to train model totally. The benchmark dataset can be defined as:
S = SIP ∪ SU (3)
where SIP represents the initial positive set containing above 58 functional AS events, and SU is the unlabeled set, containing 9,997 randomly selected ES AS events of unknown function. The SU was assembled by randomly selecting 400 AS events from each dataset, containing seven organs and three hematopoietic lineages of both human and mouse, while excluding the subset with ‘NA’ features.
Prediction workflow. The Random Forest (RF)-based algorithm was selected as the final model due to its solid performance in binary classification tasks, compared to Naive Bayes (AUC: 0.94), Logistic Regression (0.96), Support Vector Machine (SVM, 0.98), AdaBoost (0.88) and CART (0.90). Known for their robust performance across different datasets, RF algorithms can handle various feature types—categorical, Boolean, and continuous—without needing extensive feature selection. They effectively manage correlated features and require fewer hyperparameters to avoid overfitting. Moreover, RF models can reliably predict even with missing data and are adaptable to datasets lacking complete predictive feature coverage, making them applicable across different species and annotation sets. While non-linear models like RF used to be seen as less interpretable, recent advancements in understanding tree-based models have improved by explaining both the global impact of features and their influence on individual predictions.
The framework of the prediction model contained three steps: (i) the initial features were generated for the functional AS based on the dynamic, conservative and structural information; (ii) The unlabeled set was labeled by an ensemble PU learning process, including positive and negative labels. All labeled samples were divided into training and test sets; (iii) The Random Forest classifiers were trained based on training dataset. Best prediction model was selected with highest performance from 10-fold cross-validation.
Predictive features. The 19 features of the predictive model can be classified into three main categories:
a. Dynamic features of genes (six features) or AS (six features). This category depicts expression characteristics in lineage differentiation or organ development, including the expression breadth (the range value, R), the expression relationship with differentiation or development (the Spearman’s rank coefficient, ρ; the Spearman’s test p-value, Pρ), the expression specificity (the specificity index, τ) and the expression gradient (the linear regression coefficient, β; the p-value of the regression variable, Pβ).
b. Evolutionary Feature. This category captures cross-species sequence conservation, assessed by mapping vertebrate orthologs to each variant and counting the number of orthologs that are correctly aligned without gaps using CORSAIR85.
c. Structural Feature. Structural characteristics of isoforms were derived from APPRIS database76, including six features which are the structural homologs and integrity score, functionally important residues score, domain integrity score, trans-membrane helices score, signal peptide and subcellular location scores.

Positive unlabeled learning based on random forests. To alleviate the false negative problem of the unknown functional AS and to fully utilize the limited validated functional AS, positive unlabeled learning (PU learning)86 was adopted in this study and combined with random forests (PU-RF) to determine the potential functionality of AS. Moreover, cross-validation is a key approach to evaluate the performance of predictors. We employed the 10-fold cross-validation to evaluate the performance of PU learning result. The SP were randomly divided into 10 folds, with one-fold defined as the validation set Sval and the remaining nine folds used as the initial positive set SIP.
In step 1, the ‘initial negative’ training subset was randomly selected from the unlabeled subset SU with the same number of positive subsets SIP. Numerous learning processes were implemented based on T (T=100) different initial training datasets. For each base learning process, the initial training dataset was composed of two subsets:
Strain_i = SIP_i ∪ SIN_i i=1, 2, :::, T (4)
with |SIP_i| = |SIN_i| i=1, 2, :::, T (5)
where SIP is the initial positive training subset, and SIN is the initial negative training subset. The numbers of AS in SIP and SIN are the same. Because there is no validated non-functional AS, SIN is composed of the randomly selected AS from the ‘unknown’ subset SU.
The initial training dataset was used to train the classifier based on the support vector machine (SVM). And the rest of the unlabeled subset were identified with probability scores by the classifier, the subsets with scores greater than 0.8 were considered as ‘reliably positive’ SRP and less than 0.2 were considered as ‘reliably negative’ SRN.
In step 2, the ‘initial and reliably positive’ subset and the ‘initial and reliably negative’ subset above step 1 were used to train the SVM classifier again. The classifier was used for the rest of the unlabeled dataset. The reliably positive and negative samples were selected by the probability scores greater than 0.9 and less than 0.1. The step was iterated until the number of remaining unlabeled datasets (N) reached zero or no new reliable samples were generated.
SP_j+1 = SP_j + SRP_j j=0, 1, 2, :::, L (6)
when j=0, SP_1 = SIP (7)
SN_j+1 = SN_j + SRN_j j=0, 1, 2, :::, L (8)
when j=0, SN_1 = SIN (9)
Strain_ij = SP_ij ∪ SN_ij i=1, 2, :::, T j=0, 1, 2, :::, L (10)
where SP_j represents the total positive subset in the jth iteration in step 2, and SRP_j represents the reliably positive subset in the jth iteration in step 1. Similarly, SN_j represents the negative subset in the jth iteration in step 2, and SRN_j represents the negative positive subset in the jth iteration in step 1. The training set Strain_ij in the jth iteration based on the ith initial training datasets were composed of two subsets: SP_ij and SN_ij. Finally, classifiers with a recall of less than 0.85 on the validation set Sval were filtered out from the set of T classifiers. The unlabeled AS SU were definitively classified based on the average score from the retained classifiers. The final labeled dataset was split into training and test sets in an 8:2 ratio. A RF model was trained on the training set with mtry=7 and ntree=1000, and its performance was evaluated on the test set.
SLabel = Strain = T (11)
SLabel = Strain + Stest (12)

Performance evaluation of the model
The area under the receiving operator characteristic curve (AUROC)87, the associated precision-recall curves (AUPRC)88,89 and F1 score were used to evaluate the prediction quality of FAScore.
Sensitivity = TP / (TP+FN) (13)
Specificity = TN / (TN+FP) (14)
Precision = TP / (TP+FP) (15)
Recall = TP / (TP+FN) (16)
AUROC = ∫01 TPR(t) d(FPR(t)) (17)
AUPRC = ∫01 Precision(t) d(Recall(t)) (18)
F1 = 2 × Precision × Recall / (Precision + Recall) (19)
where TP, FP, TN and FN represent true positives, false positives, true negatives and false negatives, respectively.
To evaluate the ability of FAScore to identify functional AS, all AS in the test set were ranked in descending order according to the average scores. The rank index90–92 was used to measure the average rank of the functional AS in all the AS of test set. The rank index is defined as:
rank index = (1 / |Stest|) * Σ (r_b / |Stest|) (20)
where |Stest| is the number of all AS in the test set Stest, and |SP_test| is the number of AS in the positive test subset SP_test. r_b is the rank position of functional AS b in all the AS of test set. The lower value of the rank index represent the better performance of FAScore.

Features of importance
We applied the mean decrease in the Gini Coefficient45 to evaluate the importance of the feature. This coefficient indicates how well features split the data at each node in the trees, reflecting its substantial role in improving node purity and classification accuracy. Gini impurity reduction was calculated for each split in decision trees. The decrease in Gini of each feature was summed up across all the nodes and trees. The sum was divided by the total number of splits of features, producing the Mean Decrease Gini index for each feature. The mean decrease Gini index was calculated from the final model and transformed to log2 value.

Step 4: Calculation of the probability to generate the functional AS scores
By evaluating the performance of the model, we chose the best model with the highest AUROC value, AUPRC value, F1 score and the lowest rank index from 10-fold cross-validation. The probability calculated by the best model was the functional AS score (FAScore).
We applied a Gaussian Mixture Model (GMM) with three components using the R package mclust (v5.4.9) to model the distribution of FAScore. The parameters include G=3 and modelNames=“V”. GMM assumes that the data is generated from a mixture of Gaussian distributions, and is represented as:
p(x) = Σ k=1 3 π_k N(x | μ_k, Σ_k) (21)
where π_k represent the mixture proportions, and N(x, j, μ_k, Σ_k) is the Gaussian distribution with mean μ_k and covariance Σ_k for each component k.
Based on GMM, AS events were classified to three classes including functional AS, non-functional AS and uncertain AS.

Verification of functional AS. The predicted functional AS events were validated for their significance from four different perspectives, including conservation of splice site flanking sequence, function enrichment analysis of host genes, protein functional domain and post-translational modification analysis in regions of alternative splicing exons.
(1) Conservation of splice site flanking sequence
The UCSC phastCons files store conservation scores for the human genome (hg38.phastCons100way.bw) calculated from multiple alignments with other 99 vertebrate species and for the mouse genome (mm10.60way.phastCons.bw) calculated from multiple alignments with other 59 vertebrate species93. Two PhastCons files were downloaded from the UCSC database. The PhastCons conservation scores of upstream 100bp region (intronic) and downstream 50bp region (exonic) of splicing sites for the functional and non-functional AS were compared.
(2) Enrichment analysis
Gene Ontology (GO) enrichment analysis was performed separately for each lineage using the R package clusterProfiler94 (v4.2.2). For each lineage, the set of genes harboring functional and uncertain ES events was tested for enrichment against a background set comprising all genes that were expressed (TPM>=0.5) in at least one sample within one lineage to control for cell lineage-specific expression biases. And “GSEA” function was used to perform enrichment analysis of cell lineage-specific genes harboring all AS events from the differentiation.
(3) Protein functional domain and post-translational modification analysis
The protein domains and motifs of human and mouse were obtained from the UniProt database (http://www.uniprot.org/)48. The functional information was classified by Singh95, containing the annotation information from the UniProt, including active site, domain, transmembrane region, repeat, zinc finger region, compositionally biased region, DNA-binding region, region of interest, lipid moiety-binding region, short sequence motif, calcium-binding region, nucleotide phosphate-binding region, metal ion-binding site, topological domain. These regions were further classified into seven categories:
a. active site and catalytic site;
b. DNA binding domain (DBD), including C2H2-type, PHD-type, C3H1-type, KRAB, Bromo, Chromo, DNA-binding, C4-type, CHCR, A.Thook, bZIP, bHLH, CCHC-type, CHCH, Bromodomain-like, CH1, C6-type, A.Thook-like, C4H2-type and CHHC-type;
c. protein-protein interaction domain (PPI), including WD, ANK, TPR, LRR, HEAT, Sushi, EF-hand, ARM, PDZ, PH, SH3, RING-type, LIM zinc-binding, WW, SH2, BTB, FERM, CH, Rod, Coil1A, MH2, WD40-like repeat, t-SNARE coiled-coil homology, Coil1B, Cb1-PTB, Coil, CARD, SH2-like, DED, IRS-type PTB, SP-RING-type, EF-hand-like, RING-CH-type, v-SNARE coiled-coil homology, Arm domain, LIM protein-binding, GYF, PDZ domain-binding and PDZD11-binding;
d. RNA binding domain (RBD), including RRM, SAM, KH, DRBM, RBD, Piwi, PAZ, S1 motif, Punmilio and THUMP;
e. transmembrane domain (TMD), including transmembrane region, ABC transmembrane type-1, ABC transporter and ABC transmembrane type-2;
f. Signal peptide;
g. Intrinsically disordered region (IDR), disordered.
(4) Post-translational modification (PTM) analysis
The PTM information of human and mouse was obtained by the PhosphoSitePlus database49 including acetylation (AC), methylation (ME), phosphorylation (PP), sumoylation (SUMO), ubiquitination (Ubiq). The position of the above region and PTM were converted to the genome coordinate position, and intersected the position of the functional or non-functional AS.

Analysis of homologous exonic parts. Using liftOver96 from the UCSC database, the exonic part regions of human were as used input to retrieve the corresponding homologous exonic part regions (EPRs) from other species (rhesus, mouse, rat, rabbit and zebrafish) in FHO datasets. EPRs with a single corresponding region mapped in other species were retained for downstream analysis. The homologous alternative exonic part (AS exon) was defined to undergo alternative splicing in the species. Based on the evolutionary relationship of the species involved in this study, the exon evolutionary ages of human were inferred with parsimony97 For a given AS exon in human, we obtained the exon family to which this exon belong and determined the species with the greatest evolutionary distance from human among those covered by this exon. The evolutionary distance between this species and human is used as the exon evolutionary ages.

Sequence alignment and phylogenetic construction of TBC1D23. The DNA and amino acid sequences of long isoform (TBC1D23-L) were downloaded for multiple species, including Homo sapiens (ENST00000394144), Macaca mulatta (ENSMMUT00000028144), Mus musculus (ENSMUST00000023431), Bos taurus (ENSBTAT00000032386), Oryctolagus cuniculus (ENSOCUT00000017705), Gallus gallus (ENSGALT00000024650), Xenopus tropicalis (ENSXETT00000025919), Danio rerio (ENSDARG00000044357), Caenorhabditis elegans (F20D1.2.1) and Branchiostoma floridae (XP_035681897.1) from ENSEMBL or NCBI databases. Multiple sequences were aligned using Jalview (2.11.2.4) based on Mafft algorithms with default parameters. The phylogenetic tree was constructed by MEGA (v11) with Neighbor-Joining algorithms with default parameters. Based on the amino acid alignment, the sequences of the 15th exon of TBC1D23 and the flanking regions were extracted from the genomes of multiple species.

Pattern clustering of FAScores of AS. The FAScore of the ancient AS from age class 5 were extracted, and the median of all species was used to cluster functionality patterns by ClusterGVis package (v0.1.1, https://github.com/junjunlab/ClusterGVis) choosing cluster method Mfuzz. The matrix was filtered with a minimum standard deviation threshold of 0.1, and the number of clusters was inferred from the retained AS events by the getClusters function (clusters = 24). The cluster results were visualized using the visCluster function.

Mice. C57BL/6J mice (stock no. N000013) were purchased from GemPharmatech Co. Ltd. The Tbc1d23ΔE15/ΔE15 (E15-KO) mice were generated in the C57BL/6J background by GemPharmatech Co. Ltd, whose introns flanking exon 15 of Tbc1d23 were targeted by 4 sgRNAs using CRISPR/Cas9 technology. Rosa-Cas9 knock-in mice were originally purchased from the Jackson Laboratory (stock no. 024858) and backcrossed onto C57BL/6J background. Both heterozygous (Cas9+/-) and homozygous (Cas9+/+) Rosa-Cas9 knock-in mice were used in this study and collectively designated as Cas9+ (Supplementary Fig. 6d). All mice used were female and were between 8 and 12 weeks of age. The oligonucleotide (oligo) sequences of sgRNA and genotyping primers were listed in Supplementary Data 11. Unless otherwise stated, all oligos used in this study were synthesized by Sangon Biotech.
All mice were housed and bred under SPF (Specific Pathogen Free) conditions with a 12hr light/dark cycle at the Laboratory Animal Center of West China Second University Hospital, with free access to food and water. All animal experiments were carried out following the protocols approved by the ethics committee of West China Second University Hospital [(2018) Animal Ethics Approval No. 004].

Genotyping. The Mouse Tail Direct PCR kit (FOREGENE, #TP-01331) was used to extract DNA from mouse tails. Genomic DNA was extracted from sorted GFP+ cells using the TIANamp Genomic DNA Kit (TIANGEN, #DP304) and cells were lysed in elution buffer containing proteinase K and incubated at 65°C for 1 hour, followed by heat inactivation at 95°C for 5 minutes. PCR was performed for genotyping with primer sequences were listed in Supplementary Data 11.

Zebrafish. All zebrafish (Danio rerio) experiments were performed according to standard procedures as previously described56,98,99. Briefly, both adult fish and embryos were maintained at 28.5°C in Aquatic Ecosystems. AB (wild-type) and Tg(gata1a:DsRed) strains were used in this study. Larvae were imaged at 48 h post fertilization (hpf), a developmental stage at which sex is not yet distinguishable. All experimental protocols were approved by the Animal Ethical Committee, West China Hospital of Sichuan University.

Morpholino injection and validation. Antisense morpholino oligonucleotides (MO) were purchased from Gene Tools. Tbc1d23-E15-MO (E15-MO) targets the exon 15 of the zebrafish Tbc1d23 gene. Control-MO (Ctrl-MO) is a standard, mismatched control. MO were injected (3ng per embryo) into the embryo at the one-cell stage. To determine the efficiency of MO, total RNAs of ~40 zebrafish embryos at 48 hpf were extracted and were reverse transcribed to complementary DNA (cDNA) for RT-qPCR (qPCR) and RT-PCR (see below). Oligo sequences of MO and primers were listed in Supplementary Data 11.

Zebrafish fluorescent imaging. Live embryos were anesthetized and photographed. The fluorescent images were captured by a ZEISS AXIO Zoom (V16) microscope and analyzed using the ZEN (3.1) software and ImageJ (1.52a, github.com/imagej/ImageJ).

Whole Embryo Staining for Globin Expression. O-dianisidine staining was performed to determine the hemoglobin level in the zebrafish embryos as described before100,101. The dechorionated embryos at 48 h post-fertilization (hpf) were incubated in the O-dianisidine (Sigma, #243-737-5) staining solution (100 mg O-dianisidine in 70ml 100% ethanol) for 30 min in dark at room temperature (RT). Embryos were then washed with PBS with 0.1% Tween 20 (Sigma, #9005-64-5) to stop the reaction and fixed in 4% paraformaldehyde (Biosharp, #BL539A) for at least 3 h before storage in PBS at 4°C. The stained embryos were imaged on a Nikon SMZ18 stereomicroscope, and the images were used to quantify the staining signals using ImageJ.

RNA extraction, cDNA synthesis and PCR. Total RNAs of zebrafish embryos, mouse organs or cell lines were extracted using TRIzol (Life Technologies, #15596026) according to the manufacturer’s instructions. Complementary DNA (cDNA) was synthesized using the ReverAid first-strand cDNA synthesis kit (Thermo Scientific, #K1621). RT-PCR was performed using the 2× Phanta UniFi Master Mix (Vazyme, P516). And real-time qPCR was performed using the Blue qPCR SYBR Master Mix (YEASEN, #11184) on a Bio-Rad (CFX Connect) thermocycler. Primers were listed in Supplementary Data 11.

Blood routine, cell isolation and flow cytometry analysis. The blood routine was assayed using BOWLINMAN (BM660). Whole BM cells were collected from the tibia, femur and ilium, filtered through 70-μm nylon mesh, quantified by PBS and counted by Countstar Altair (Counstar). For the analysis of hematopoietic stem and progenitor cells (HSPCs), BM cells were lysed with red blood cell lysis buffer (Solarbio, #R1010), washed with PBS and filtered through 70-μm nylon mesh.
HSPCs were further enriched using c-kit/CD117 MicroBeads (Miltenyi Biotec, #130-09224) or EasySepTM Mouse Hematopoietic Cell Isolation Kit (STEMCELL Technologies, #19856) with a biotinylated anti-mouse lineage cocktail (CD11b, B220, Gr-1, TER119 and CD3e) according to the manufacturer’s instructions.
Single cell suspensions were stained with antibodies at 4°C for 30 min in dark. The antibodies used in the study were listed in Supplementary Data 12. Flow cytometric analysis was performed using an Attune Nxt (Life Technologies) flow cytometer and cell sorting were performed using a FACSAria III (BD Biosciences). The data were analyzed using FlowJo software (v.10.8.1, BD). Gating strategies for all flow cytometry experiments are illustrated in Supplementary Fig. 13.

Colony-forming unit assay. Freshly isolated whole BM cells, enriched c-kit+ cells or sorted GFP+ cells were plated in methylcellulose culture medium (MethoCult™, STEMCELL Technologies) and incubated at 37°C with 5% CO2. For cells plated on MethoCult™ GFM3334, colonies were counted and scored as colony-forming unit erythroid progenitor cells (CFU-E) at 48h after plating. Cells plated on MethoCult™ GF M3436 were evaluated 10–14 days after plating and scored as burst-forming unit erythroid progenitor cells (BFU-E). For cells plated on MethoCult™ GFM3434, colonies derived from BFU-E and granulocyte-macrophage progenitor cells (CFU-GM, CFU-G and CFU-M), multipotent granulocyte, erythroid, macrophage, megakaryocyte progenitor cells (CFU-GEMM) were counted 9–12 days after plating.
For cell isolation from MethoCult™ GF M3434 colonies, the methylcellulose medium containing colonies was collected. The sample was washed five times with 10 volumes of PBS. After centrifugation, the supernatant was removed to eliminate residual methylcellulose. Cells were stained with antibodies and sorted by flow cytometry. Gating strategy was illustrated in Supplementary Fig. 13a.

Single-cell mRNA sequencing. Female E15-KO and the WT mice aged 8–9 weeks were used for enrichment of LIN- HSPCs as mentioned above. The collected cells were resuspended in PBS containing 1% BSA at the concentration of 1×106 cells/mL and prepared for library preparation. Single cells were prepared in the Chromium Single Cell Gene Expression Solution using the Chromium Single Cell 3’ Gel Bead, Chip and Library Kits v2 (10× Genomics) according to the manufacturer’s protocol. Briefly, the cells were loaded in each channel and then partitioned into nanoliter-scale Gel Bead-In-EMulsions (GEMs) using the Chromium instrument, where all cDNA generated from an individual cell share a common 10× Barcode. To identify the PCR duplicates, the Unique Molecular Identifier (UMI) was also added. The GEMs were incubated with enzymes to synthesize full-length cDNA, which was subsequently amplified by PCR to generate a sufficient quantity for library construction. After cDNA amplification, enzymatic fragmentation and size selection were performed using SPRI select reagent (Beckman Coulter, #B23317) to optimize the cDNA size. P5, P7, a sample index and read 2 (R2) primer sequence were added by end repair, A-tailing, adaptor ligation and sample-index PCR. Amplified cDNA and final libraries were assessed with the Agilent BioAnalyzer using a High Sensitivity DNA Kit (Agilent Technologies). Finally, the libraries were sequenced on the Illumina NovaSeq 6000 platform by Novogene, Beijing, China. On average, 101 Gb of raw data were generated for E15-KO and 88 Gb for WT HSPCs.
The sequencing data were processed using the Cell Ranger (v6.1.2) with default parameters and the GRCm38 mouse reference genome (downloaded from https://ftp.ensembl.org/pub/release-93/fasta/mus_musculus/dna/). The raw gene expression matrices were merged using R (v4.0.2) and converted to a Seurat object using Seurat (v4.0.2). The doublets were identified using DoubletFinder (v2.0.3). Cells having less than 200 unique molecular identifiers (UMIs) or more than 20% mitochondrial counts, or genes below 500 UMIs in all cells were filtered. The batch effect was corrected by a single integrated analysis.
Gene expression matrices of the 12,720 cells after the filter were normalized to the total cellular read count using linear regression using the ScaleData function in Seurat. To reduce the dimensionality of this dataset, the top 2000 highly variable genes were summarized using principal component analysis (PCA). Uniform manifold approximation and projection (UMAP) was generated from the first 20 dimensions using RunUMAP with default parameters. Cell types or clusters were annotated according to the canonical cell markers.
The differentially expressed genes (DEGs) between E15-KO and WT in each cell type were identified using the FindMarkers function in Seurat with default parameters (log2 fold change ≥ 1, and Bonferroni-adjusted p value < 0.05, Wilcoxon Rank Sum test). Cellular differentiation states were predicted by using CytoTRACE (v0.3.3). The Augur (v1.0.3) was used to identify which lineage exhibit the most transcriptional changes between E15-KO and WT. Results were visualized using the ggplot2 package (v3.3.0).
Co-expression analysis was carried out using the hdWGCNA R package60 (v0.2.2). Erythroid-related cells including ProE/BasoE, PolyE, Early-orthoE and Late-orthoE were selected for co-expression network analysis. We retained 19,655 genes which were expressed in at least 5% of cells from any cluster. Metacell transcriptomic profiles were constructed separately for each of the 2 conditions and each cell type using the hdWGCNA function MetacellsByGroups with parameters ‘k=10, max_shared = 10, min_cells = 0’. We selected a soft-power threshold β=16 based on the parameters we performed with the TestSoftPowers function. The co-expression network was computed with the ConstructNetwork function with the default parameters. Module eigengenes were computed using the ModuleEigengenes function, and we applied Harmony to correct MEs based on the sequencing batch. Eigengene-based connectivity for each gene was computed using ModuleConnectivity. The co-expression network was embedded in two dimensions using UMAP with the RunModuleUMAP function with the top five genes (ranked by eigengene-based connectivity) as the input features. Differential module eigengene (DME) analysis was performed using the FindDMEs function, and the Wilcoxon test was used to compare KO and WT groups.
The co-regulatory networks of transcription factors and potential target genes were inferred using pySCENIC62 (v0.12.1) from scRNA-seq. The erythroid-related cell types were retained. The genes with a total count of less than 27 (1 UMI x 1% of cells) or expressed in less than 1% of the cells are filtered out. 11,450 genes were retained to construct of gene regulatory network using ‘pyscenic grn’ with default parameters. And we chose the motif database file with ranking in which matrix containing motifs as rows and genes as columns and ranking position for each gene and motif (based on CRM scores) as values. The search space around the TSS of the gene in which the motif is scored is defined as 500 bp upstream of the TSS and 100 bp downstream. The motif annotations file of mice based on the 2017 cisTarget motif collection was downloaded. We used the ‘pyscenic ctx’ function to predict the regulons and find the enriched motifs with parameters ‘--mode dask_multiprocessing’ and ‘--mask_dropouts’. In the following, the activity of regulons was quantified by calculating AUC using the ‘pyscenic aucell’ function. The co-regulatory network of TFs and their target genes was visualized using Gephi (v0.10.1)102.

Cell culture and transient transfection. HEK-293T and K562 was purchased from Type Culture Collection of the Chinese Academy of Sciences. HEK293T cells were cultured in DMEM medium (Hyclone, #SH30022.01) and K562 cells were cultured in RPMI 1640 medium (Hyclone, #SH30809.01), both supplemented with 10% fetal bovine serum (FBS, NEWZERUM, #FBS-S500) and 1% penicillin-streptomycin (Hyclone, #SV30010).
c-kit+ cells of WT or Cas9+ mice were enriched and cultured in 50% IMDM (Gibico, #12440053) and 50% DMEM medium supplemented with 20% FBS (Gibico, #10099141C), 10 ng/mL stem cell factor (SCF, Peprotech, #250-03-10), 10 ng/mL interleukin-6 (IL-6, Peprotech, #216-16-10), 2 ng/mL interleukin-3 (IL-3, Peprotech, #213-13-10), 1% penicillin-streptomycin, and 50μM beta-mercaptoethanol (Sigma, #SLM–063). All cells were incubated at 37°C in 5% CO2. After 14–16 h, the cultured cells were transduced with retrovirus particles.
Transient transfections were performed with overexpression plasmids (described below) using the Liposomal Transfection Reagent (YEASEN, #40802ES03). The medium was replaced at 4–6h after transfection, and the cells were harvested at 48h after transfection.
For viral transduction, HEK293T cells were plated onto six-well plates in 2ml complete DMEM medium. The next day, cells were transfected with a confluence of 80–90%. Lentiviruses particles were generated by CaPO4-mediated transfection of HEK293T cells with 2µg of psPAX2 (Addgene, plasmid #12260), 1μg of pMD2.G (Addgene, plasmid #12259) packaging plasmid, and 4μg of the lentiviral vector plasmid. Alternatively, lentiviral particles were generated by Liposomal Transfection Reagent transfection of HEK293T cells with 1µg pMD2.G, 1µg pMDLg/RRE (Addgene, plasmid #12251), 1µg pRsv-Rev (Addgene, plasmid #12253) packaging plasmid and 3µg of the lentiviral vector plasmid. Retroviruses particles were generated by CaPO4-mediated transfection of HEK293T cells with 1.5µg pCL-Eco (Addgene, plasmid #12371), 0.2µg pVSV-G (Addgene, plasmid #138479) packaging plasmids and 4µg Retroviral vector plasmid. The medium was exchanged for fresh complete medium 6h after the transfection. All viral supernatants were harvested at 48h and passed through a 0.45μm filter to remove cells.
HEK293T cells were transduced with lentivirus particles in the presence of 4μg/mL polybrene (YEASEN, #40804ES76). K562 cells were transduced with lentivirus particles and spun at 1000g at 37°C for 1h in the presence of 4μg/mL polybrene. c-kit+ cells were transduced with retrovirus particles and spun at 1,000g at 37°C for 1h in the presence of 8μg/mL polybrene.

Generation of cell lines expressing only TBC1D23-S or TBC1D23-L isoform. We first generated a TBC1D23 knock-out (KO) line by CRISPR/Cas9 methodology, using lentiCRISPRv2 (Addgene, plasmid #52961). The resulting constructs were co-transfected with lentiviruses packaging plasmids (psPAX2 and pMD2.G) to generate lentiviruses particles as described above. The viral supernatants were harvested and added into K562 or HEK293T cells with 4μg/mL polybrene. K562 or HEK293T cells were selected with 2μg/mL or 1μg/mL puromycin (BBI lifesciences, #A600013), respectively at 72h after infection. Cloning was performed using limiting dilution in 96-well plates and the disruption of TBC1D23 was confirmed using Sanger sequencing of the targeted region and immunoblotting using anti-TBC1D23 antibodies.
For consecutive expression, TBC1D23-S and TBC1D23-L were amplified from K562 cDNAs and cloned into pLVX-mcherry-N1 (TaKaRa, #632562) by Easy Clone Kit Ⅰ(YoungGen, #CL101A), resulting in overexpression constructs having N-terminally mCherry-tagged TBC1D23-S or TBC1D23-L under a CMV promoter. The overexpression constructs were introduced into the cloned TBC1D23 KO line by lentiviral transduction as described above and the cells were selected with 1mg/mL G418 (BBI lifesciences, #A600013). The expression TBC1D23-S and TBC1D23-L was confirmed using immunoblotting (Supplementary Fig. 9a). Oligo sequences for gRNA or primers were listed in Supplementary Data 11.
For inducible expression, TBC1D23-S and TBC1D23-L were cloned into the tetracycline-inducible lentiviral expression vector pTSB-Tight-EF1-tetR-F2A-Puro, which is a tetracycline (Tet)-inducible expression system103. While this system enables high-level transgene expression via Tet operator binding in the presence of tetracycline, the mini-CMV promoter lacks enhancer elements and consequently drives only low-level basal expression in the absence of tetracycline. Lentiviral particles were produced as above and used to transduce TBC1D23-KO HEK293T cells. The expression levels of TBC1D23-S and TBC1D23-L were modulated by titrating the doxycycline concentration to have the comparable TBC1D23 protein level as the control cell line. The oligo sequences for gRNA or primers were listed in Supplementary Data 11.

Generation of overexpression cell lines. SUMO1 was amplified from HEK293T cDNAs and cloned into pLVX. The resulting overexpression constructs contain an N-terminally HA-tagged SUMO1 expression cassette and neomycin selection marker. The overexpression constructs were introduced into HEK293T cells by lentiviral transduction with lentiviral packaging plasmids (pMD2.G, pMDLg/RRE, pRsv-Rev) as described above. The target cells were selected with 500μg/mL G418. The cDNA fragments of HDAC1 and RanGAP1 were amplified from the cDNAs of K562 cells. HDAC1 was cloned into pcDNA3.1(+) with an N-terminal FLAG tag. Mutant HDAC1 (HDAC1mut, K444R/K476R) was generated by PCR mutagenesis using 2× Phanta UniFi Master Mix. RanGAP1 was cloned into pcDNA3.1(+) with an N-terminal HA tag. A DNA fragment containing the N-terminally GFP-tagged E3 ligase domain of RanBP2 was synthesized (GENEWIZ) and cloned into PcDNA3.1(+). The E3 catalytic dead mutation of RanBP2 (F2735A/F2736A/C2737A) was generated by PCR mutagenesis using the High-Fidelity PCR kit (MACLAB, #I5HM). The primers were listed in Supplementary Data 11. These overexpression constructs were transiently introduced into HEK293T cells or HA-SUMO1-expression HEK293T cells.

Generation of exon-deleted and overexpression mouse HSPCs. The target exon of Tbc1d23, Ebp41l1, Ssbp3, and Klf6 were deleted by two sgRNAs (listed in Supplementary Data 11) were cloned into pMSCV-U6-sgRNA-IRES-GFP, which was modified from pMSCV-IRES-GFP (Addgene, plasmid #20672). The coding sequence of HDAC1 and HDAC1mut were cloned into pMSCV-IRES-GFP for overexpression. The resulting constructs were cotransfected with retroviral packaging plasmids (pCL-Eco, pVSV-G) to generate retroviral particles as described above. The viral supernatants were harvested and added into c-kit+ cells of Cas9+ or WT mice with 8μg/mL polybrene. At 24h after transfection, the target cells were enriched by GFP.

Immunoprecipitation and Immunoblotting. Cells were harvested and lysed in RIPA (Thermo Fischer Scientific, #89900) containing Mini EDTA-free protease inhibitor cocktail (Thermo Fischer Scientific, #A32955) on ice for 30min, and pelleted. The supernatant was solubilized in 0.1 M Tris-HCl containing 2% SDS and 0.05 M DTT at 95% for 8min.
For immunoprecipitation (IP), cells were harvested and lysed in Western-IP lysis buffer (Beyotime, #P0013) containing 10mM protease inhibitors (Bimake, #B14002), 20mM N-Ethylmaleimide (Sigma, #128-53-0), 1mM DTT similar to previous studies104. Lysates were centrifuged at 14,000g at 4°C for 20min. Rabbit anti-TBC1D23 (Proteintech, #17002-1-AP), mouse anti-FLAG (Proteintech, #20543-1-AP), Rabbit anti-HDCA1 (CST, #34589) antibodies were equilibrated with Tris-buffered saline (TBS) (Sangon Biotech, #C520002) before use. The supernatant from cell lysates was mixed with the protein A/G magnetic beads (Bimake, #B23202) and incubated with antibodies at 4°C overnight. The protein-antibody-magnetic beads complexes were then washed using TBS three times and eluted with sodium dodecyl sulfate (SDS) loading buffer (YEASEN, #20315ES) at 95°C for 10min. The solutes were loaded and separated by 10% SDS–polyacrylamide gel electrophoresis (SDS-PAGE) gel (Epizyme, #PG112) for mass spectrometry or immunoblotting.
Protein extracts or eluted proteins from the IP assay were separated on 10% SDS polyacrylamide gels and transferred to 0.45μm polyvinylidene difluoride (PVDF) membranes (Beyotime, #FFP28). The membranes were blocked in TBS buffer with 0.1% Tween 20 and 5% skimmed milk at RT for 1hr. The blots were probed with Rabbit anti-TBC1D23 (1:1000), Mouse anti-HDCA1 (1:1000) (CST, #5356), Rabbit anti-sumo1 (1:1000) (HUABIO, #ET1606-53), Rabbit anti-Sumo2/3 (1:1000) (HUABIO, #ET1701-17), Mouse anti-Flag (1:1000), Mouse anti-GFP (1:1000), Mouse anti-HA (1:1000) or Mouse anti-RanBP2 (1:1000) (Santa Cruz Biotechnology, #sc-74518) at 4°C overnight. Then the membranes were incubated with the Goat anti-Rabit IgG Seconadary antibody HRP conjugated IgG (1:2000) (SAB, #L3012-2) or Goat anti-Mouse IgG Seconadary antibody HRP conjugated IgG (H+L) (1:2000) (SAB, #L3032-2) at RT for 1h. The chemiluminescence signal was detected using a ChemiDOC™ MP Imaging System (Bio-Rad Laboratories). ImageJ software was used to analyze the gray value of the WB.

IP-Mass. Gel pieces (number of technical and/or biological replicates=1/group) were destained in 50mM NH₄HCO₃/50% ACN, dehydrated with ACN, and subjected to reduction with 10mM DTT (56°C, 60min) and alkylation with 55mM iodoacetamide (room temperature, dark, 45min). After washing, gel pieces were digested with 10ng/μL trypsin at 37°C overnight. Peptides were extracted with 50% ACN/5% FA followed by 100% ACN, dried, and reconstituted in 2% ACN/0.1% FA.
The tryptic peptides were dissolved in 0.1% formic acid (solvent A), directly loaded onto a home-made reversed-phase analytical column (15-cm length, 75μm i.d.). The gradient was comprised of an increase from 6% to 23% solvent B (0.1% formic acid in 98% acetonitrile) over 16min, 23% to 35% in 8min and climbing to 80% in 3min then holding at 80% for the last 3min, all at a constant flow rate of 400nl/min on an EASY-nLC 1000 UPLC system. MS analysis was performed on a Q Exactive™ Plus mass spectrometer with NSI ionization (2.0kV). Full MS scans (m/z 350–1800) were acquired at 70,000 resolutions, followed by data-dependent MS/MS scans (NCE 28, 17,500 resolution) with a dynamic exclusion of 15.0s. Data were processed using Proteome Discoverer search engine (v. 2.4). Tandem mass spectra were searched against Homo Sapiens (SwissProt, 20366 entries) database with trypsin/P digestion (max 2 missed cleavages), precursor mass tolerance of 10ppm, fragment tolerance of 0.02 Da, carbamidomethyl (C, fixed), and oxidation (M, variable). Peptides were filtered at high confidence with ion score > 20.

Immunofluorescence. The cells were loaded on poly-lysine treated slides and fixed with 4% paraformaldehyde at RT for 15min. After washing with PBS for three times, cells were permeabilized by 0.1% Triton X-100 at RT for 15min followed by blocking with 10% FBS in PBS at RT for 1h. Cell staining was performed using primary antibodies including Rabbit anti-TGN46 (1:300, Abcam, ab50595), Mouse anti-GM130 (1:200, BD, #610822), Mouse anti-Cl-MPR (1:200, Bio-Rad, #MEM-2), Rabbit anti-ZFPL1 (1:200, Invitrogen, #PA5-53593), Rabbit anti-TBC1D23 (1:50, Sigma, #HPA038152) or Mouse anti-RanBP2 (1:1000) at 4°C overnight. The cells were washed again with PBS for three times and incubated with secondary antibodies, including FITC affinipure goat anti-Mouse IgG antibody (1:2000, Jackson ImmunoResearch), FITC affinipure goat anti-Rabbit IgG antibody (1:2000, Jackson ImmunoResearch, #115-095-003) or Alexa Fluor 647 affinipure goat anti- Mouse IgG antibody (1:2000, Jackson ImmunoResearch, #115-605-003) at 37°C for 1h. The nuclei were stained with DAPI (Sangon Biotech, #E607303). After mounting with VectaShield medium (Vector Laboratories), images were acquired using confocal microscopy (Olympus V3000). The images were used to quantify the Fluorescence intensity using ImageJ.

SUMO proteome analysis. TBC1D23 KO K562 cells, or cells expressing only TBC1D23-L or -S (cell number = 4×107, number of technical and/or biological replicates = 1) were collected for quantification of SUMO modifications using a PTMScan HS Ubiquitin/SUMO remnant kit (Cell Signaling Technologies, #59322) according to the manufacturer’s instructions. In brief, the cells were harvested in 9M Urea Lysis buffer. The lysate was sonicated using a sonicator (BioSafer, #650-92) with a microtip at 15W output with 3 bursts of 15s each. Clear the lysate by centrifugation at 20,000g for 15min at room temperature and transfer the protein extract (supernatant) into a new tube. Protein concentrations were determined by bicinchoninic acid (BCA) assay (Thermo Fisher Scientific, #23227). DTT (1.25M) was used for protein reduction, and then iodoacetamide solution (Thermo Fisher Scientific, #A39271) was added to the cleared cell supernatant, and diluted 4-fold with 20mM HEPES pH 8.0 to a final concentration of ~2M urea, 20mM HEPES, pH 8.0. The samples were digested overnight at 37°C with WALP enzyme (Cell Signaling Technology, #33036) 0.4 mg/mL stock at the ratio of 1mg enzyme: 100mg substrate. Complete digestion was confirmed by SDS-PAGE. The digested sample was acidified with 1/20 volume of 20% trifluoroacetic acid to a final concentration of 1%. The pH was checked by spotting a small amount of peptide sample on a pH strip (the pH should be under 3). After acidification, allow the precipitate to form by letting the sample stand for 15min on ice. Centrifuge the acidified peptide solution for 15min at 1,780g at room temperature to remove any precipitate. The peptide-containing supernatant was transferred into a new 50mL conical tube without dislodging the precipitated material. Acidified and cleared digest was loaded onto the C18 column, which was pre-wetted with 5mL 100% ACN and washed sequentially with 1, 3, and 6 mL of Peptide Purification Equilibration Solution (0.1% TFA). The column was then washed sequentially with 1, 3, and 6 mL of Peptide Purification Equilibration Solution (0.1% TFA) and with 2 mL of Peptide Purification Wash buffer. Place columns above new 50mL polypropylene tubes to collect eluate. Elute peptides with a sequential wash of 3mL and then 7mL of Peptide Purification Elution Solution (0.1% TFA, 50% acetonitrile). TFA was removed by freezing the eluate on dry ice (or -80°C freezer) overnight and lyophilizing frozen peptide solution for a minimum of 2 days. The tube containing lyophilized peptides was centrifuged for 5 minutes at 2,000g at room temperature to collect material for dilution in HSIAP Bind buffer #1. PTMScan HS antibody-bead slurry in ice-cold phosphate-buffered saline (PBS) was added to each clarified sample. Samples were incubated on an end-over-end rotator at 4°C for 2h. The beads were collected using a magnetic stand and then washed sequentially with HS IAP Wash buffer and LCMS water to wash the beads. The beads were resuspended in IAP Elution buffer (0.15% TFA) and the eluted sample was collected in a new microcentrifuge tube.

Mass spectrometry. Peptides were desalted and resuspended in buffer A (2% ACN, 0.1% FA). Then analyzed by EASY-nanoLC 1000 coupled to a high-resolution mass spectrometer (Q Exactive Plus, Thermo Fisher Scientific). The samples were loaded onto a 75μm (inner diameter) × 2cm (length) trap column and a 75μm (inner diameter) × 15cm (length) analytical column, which were packed with C18 particles (DIKMA) in-house. The samples were loaded into a 75μm (ID) × 2cm (length) trap column and a 75μm (ID) × 15cm (length) analytical column, respectively. Data-dependent acquisition (DDA) was performed in positive ion mode, and samples were analyzed with a 105min or 65min gradient from 12 or 13 to 100% buffer B (80% acetonitrile, 0.1% formic acid) for proteomics at a flow rate of 330nL/min. Full MS scans were obtained with a range from 350 to 1600 m/z with a resolution of 70,000 at m/z = 200, with an AGC target value of 3×106 and a maximum injection time of 20ms. The AGC value for MS/MS was set to 1×105, the maximum injection time was 100ms, and the resolution was 35,000. Precursor ions with charge state z = 1, 7–8 or unassigned were excluded. For matching and database searching methods for WaLP digested samples refer to previous studies105. GO enrichment analyses were performed by R package clusterProfiler94 (v4.2.2) for the down-regulated expressed SUMO-modified proteins of TBC1D23-KO compared to cells expressing only TBC1D23-S or TBC1D23-L.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Hu et al report on the functional landscape of alternative splicing in hematopoietic lineages. Alternative splicing is an important post-transcriptional RNA-processing step and contributes to transcript and protein isoform diversity. Importantly, while cataloging AS has become common practice determining the function impact of AS remains difficult to experimentally define and prioritize. This manuscript intends to address this gap by developing a machine learning model, FAScore, that incorporates static information (such as genomic sequence features) and dynamic features (such as expression) to predict which splicing events are functionally relevant. The model employs a random forest framework employing a list of features that incorporates expression, structure and evolutionary conservation. Using the model, they were able to characterize events into “functional” and “non-functional” categories and validate some of the “functional” ones. While there is some excitement regarding the modeling, the curation and assessment of splicing was unclear. However, the experimental follow-up in its current form is challenging to follow as it lacks details. The paper follows primarily one splicing event in the Tbc1d23-E15 in mouse\nIn summary this is a paper with a significant amount of data however in its current form is challenging to evaluate.\nMajor\nThe paper is dense and written as a short form article, but this makes it very challenging to understand and would benefit from more main figures and a significant improvement in writing (primarily adding more information).\nExperimental sections lack of sufficient information regarding species, model systems used, validation etc. Much of this information is crucial for an appropriate review.\nFor validation of mechanism, omic-scale experiments are offered such as scRNAseq, IP-mass spec, mass spec for sumo but these are provided with little context or interpretation.\nWhile I appreciate the work for obtaining FAScores, I’m not sure that the experimental components and the selection of AS for validation necessarily requires this score. One might prioritize the same events based on conservation, protein domain impact, and lineage. Thus, more should be done to discuss why the model facilitates discovery beyond the standard practices of looking at dynamics and conservation.\nDescription of “Dynamically expressed genes” is difficulty to follow. This section needs rewriting and perhaps schematics to help explain the method in more detail. Representative example data for certain genes or alternative splicing events would be helpful.\nFor expression breadth, it is mentioned that TPMs for each gene was used. TPMs are normalized within samples and are not the best measure for between-sample comparisons. Scaling it between 0and 1 could lead to confounding effects. What does the overall distribution of dynamic scores look like for genes and AS events. How was 0.4 chosen?\nFor expression gradient, it’s important to show that the relationship was linear between time points and gene expression and AS events which might not always be the case. How were these dealt with? Furthermore, no reasoning was provided for capping a continuous scale such as regression coefficient to -1 and 1. Examples of events that agree with expression specificity score would be helpful.\nAlternative splicing patterns are already shown to be lineage specific (Merkin 2013 Science), so dynamic AS being lineage specific could just be a sub-effect of that.\nSequence conservation is part of the formula used to calculate the score and is one of the top features based on gini index. Could the model be biased towards events with high conservation (i.e. is this circular)? Conservation and protein domain composition is a key feature of gene importance, thus is an ML approach required for identifying these events and prioritizing them?\nDoes any non-functional or uncertain event have literature on their functional relevance? What is the false negatives on this model?\nThe statement: \"In short, alternative splicing of TBC1D23-E15 leads to differential binding affinity to the RanBP2/RanGAP1 complex and may affect cellular SUMOylation level.\" is not fully supported by experiment. The raw gels indicate minimal difference in interaction and exclusively in an over-expression setting, further the term affinity is used but not measured using biophysical assays. Similarly, the impact of sumoylation needs to be further validated beyond Co-Ip experiments.\nMethods are lacking as is model validation for zebrafish and mouse.", "reviewer_2": "This manuscript describes a Random Forest based approach to score alternative splicing (AS) events based on their likelihood to have functional relevance for specific hematopoietic lineages. Based on this score (FAScore) as well as the differential regulation of sequence inclusion, they then selected four AS events, which they functionally validated. They follow up on one of them (exon 15 of TBC1D23), and elucidated the mechanistic details on how the skipping/inclusion of that exon leads to its physiological impact.\nOverall, I found the experimental tests and validations to be very well executed, specially for TBC1D23. The development of the FAScore is not particularly original, perhaps with the exception of making it lineage-specific (instead of \"global\"). However, the authors seem to then average FAScores and used them in combination with the DyScore (i.e. classic delta PSI) to select the events of interest (Fig. 2d).\nI have some issues/concerns:\n1) The positive training set to develop the FAScore is of 63 events according to Sup Table 3 (even though the Methods section says 71). This number seems quite low for an RF approach. Also, looking at the references in Table S3, for some of these events I am not sure there is actual evidence of functionality, only of protein translation by Western blot. This should be clarified. Finally, most of the events used for training are of the exon skipping type (see point 2 below). These limitations should be explicitly acknowledged in the main text.\n2) It is unclear to me how the benchmarking results in Figures 1 and 2 and in the associated Supplementary Figures look for each type of AS event separately. All events seem to be analyzed and treated together, which I believe is potentially problematic. As mentioned above, FAScore has been mainly trained for exon skipping events and it is known that the level of conservation, functionality, etc. is very different for different types of events.\n3) Another aspect I found unclear is how the authors have dealt with constitutive, cryptic and \"real\" AS exons. Based on the number of AS events analyzed, it seems to me they have applied their score to all exons in the genome. What does it mean if a constitutive exon has a high FAScore? Does it mean that its mis-splicing will have a functional impact (as expected)? Again, it would be good to see how the FAScore ranges for the different types of events (constitutive, AS, etc.), and how the benchmarkings perform separately per inclusion type. For instance, if all highly regulated exons (high DyScore) have high FAScores, this means the FAScore is mainly a non-useful score on its own. This score seems to categorize thousands of exons/events as functional, which is likely questionable.\n3) Perhaps more important: the standard in the field of AS to rank candidates for experimental validations is to select evolutionarily conserved events that change their inclusion levels across lineages and/or a time course. Is the FAScore approach really much better than this? The authors selected four events with high FASscores and high DyScores (i.e. change in their inclusion levels across lineages and/or the time course). If I understood correctly, the most important contributor to the FAScore is conservation and then lineage regulation (which makes it redundant with DyScores). In fact, based on Fig. 2d, the events selected for follow up have a fairly modest FAScore (0.55-0.70) and there are at least dozens of events with much higher FAScores. Why have those events not been prioritized? In other words, the FAScore is perhaps useful, but it seem potentially an oversell.\n4) What gene sets were used as background for the GO enrichment analyses? In the Methods section it is unclear if they have used any background at all. If not, it is possible that the enrichments they obtained were due to lineage-specific gene expression biases, which are known to be a major confounder for AS analyses and could explain why they fit the lineages functions well.\n5) The zebrafish morpholino experiment is much weaker than all the other experiments. In my opinion, it is OK to keep these results, but I would move them after the first results in mouse (i.e. \"white blood cell counts remained unchanged\") and introduce them in a very specific manner, e.g.: \"To assess the conservation of this phenotype, we ....\"\nIn summary, I believe the experimental and mechanistic aspects of this manuscript are much stronger than the FAScore, which is perhaps not needed. My recommendation is to downplay this part, and simply sell it as a way to prioritize exons for experimental follow up. As it is now, it has potential flaws and its relevance is oversold, which may require substantial work to solve.", "reviewer_3": "Hu et al. surveyed the transcriptome of the hematopoietic system in six vertebrate species and trained a machine learning model to identify functional lineage-specific alternative splicing events. Their approach relies on a Random Forest model trained on a large dataset of differentially spliced events between lineages and leveraging positive unlabeled learning to take advantage of a set of 71 validated functional events curated from the literature. The interpretation of the model yielded a set of 19 features which were predictive for alternative event functionality. The model was then applied to identify several candidate functional events in lymphoid, myeloid and erythroid commitment. The authors go on to show that one of their candidates, a vertebrate cassette inclusion event in TBC1D23, is important for erythropoiesis using in vivo and ex vivo experiments in the zebrafish and mouse hematopoietic cells. They finally provide mechanistic evidence to demonstrate that the inclusion of the alternative TBC1D23 exon 15 in the erythroid lineage decreases binding to its partner complex RANBP2/RANGAP1, likely leading the activation of transcription factors involved in erythroid differentiation through lower SUMOylation of HDAC1.\nThis study tackles the extent and the role of alternative splicing during hematopoiesis. This is an important and poorly studied problem with wide ranging applications in fundamental and translational immunology. Notably, the manuscript presents three valuable components to the field: 1) a detailed survey of alternative splicing events across hematopoietic lineages, between species and between developmental stages; 2) a predictive model for splicing event functionality; 3) the characterization of a novel alternative event in TBC1D23 that is key during erythropoiesis.\nThe methodology leading to the development of the FAScore and DyScore is reasonable and robust. The methods for the experimental validation of candidate functional events are generally sound, but key experiments and details lack in some sections, in particular in demonstrating the lineage specific of TBC1D23 exon 15 inclusion and differential SUMOylation of HDAC1, as described below.\nSpecific comments:\nThe astronomical number of differential/lineage-specific splicing events reported and differentially expressed genes in the TBC1D23-E15KO mouse is hard to believe and likely obscures some of the most biologically relevant events. This could be due to lenient filtering criteria and/or specific tools and methods used for transcriptomic analyses. Depending on the version and parameters used, SUPPA2 can have a very high false-positive rate. The choice of performing scRNA-Seq rather than FACS followed by bulk RNA-Seq to identify DEGs in the TBC1D23-E15KO is puzzling. The authors should at least comment on these points and possibly consider revisiting their filtering criteria.\nA lot of the evidence cited to support assumptions is anecdotal. For example, the identification of FYN-E5 and PIK3R1-E1 among high FAScore and DyScore events is used as a demonstration that high FAScore and DyScore can be used to find lineage-specific functional events, but how many known lineage-specific functional events can be found among low FAScore/DyScore events?\nThe experiments in Figure 2E-G must be better described. Currently it is impossible to tell from the main text or figure legend that these experiments use mouse progenitor cells that were transduced with exon targeting gRNAs and then sorted by lineage to perform the colony forming assays. Even the methods do not specify that the experiment was performed in mouse cells and what mouse strain were used. There should also be a panel showing baseline exon inclusion levels by RT-PCR in the control GM, GEMM and BFU-E lineages.\nThe specificity of TBC1D23 exon 15 inclusion is not shown in any of the figures. This is a critical piece of missing data. Evidence should be provided in the form of RT-PCR on several sorted human and/or mouse hematopoietic cell populations.\nDoes the MO in Fig 3A-B also affect the levels of the TBC1D23 exon skipping transcript? Currently, it is impossible to confirm that the phenotype is due to loss of TBC1D23-E15 specifically or loss/decrease of all TCB1D23 expression altogether.\nSupplementary Fig 5D should include a panel showing TBC1D23 isoform ratios by RT-PCR using a single primer pair in flanking exons.\nNowhere is the difference in HDAC1 SUMOylation between TBC1D23-E15KO and WT shown. This is another critical piece of missing data.\nTBC1D23 E15 may alter interaction with partners, but the changes shown in co-IP experiments with RANBP2/RANGAP1 are rather small and downstream changes in HDAC1 SUMOylation are not quantified; how can the loss of a TBC1D23 isoform be expected to lead to the differential expression of close to 6,000 genes?\nLine 209: “classed” should be “classified”."}, "author_response": "Response to Reviewers\nMS Submission ID: NCOMMS-25-31345\nMS title: The functional landscape of alternative splicing in hematopoietic\nlineage commitment\nDear editors and reviewers:\nWe thank the editor and the reviewers for their valuable comments and\nsuggestions. In response to the comments, we performed additional analysis\nand experiments, revised the manuscript (the changes are shown in red font)\nand updated the accompanying figures and tables. We addressed the concerns\npoint-to-point under each comment. And the changes that we introduced into\nthe revised manuscript (MS) are shown in italics with page and line numbers in\nthe revised MS indicated.\nYours sincerely,\nLu Chen, PhD\nReviewer #1:\nHu et al report on the functional landscape of alternative splicing in\nhematopoietic lineages. Alternative splicing is an important post-transcriptional\nRNA-processing step and contributes to transcript and protein isoform diversity.\nImportantly, while cataloging AS has become common practice determining the\nfunction impact of AS remains difficult to experimentally define and prioritize.\nThis manuscript intends to address this gap by developing a machine-learning\nmodel, FAScore, that incorporates static information (such as genomic\nsequence features) and dynamic features (such as expression) to predict which\nsplicing events are functionally relevant. The model employs a random forest\nframework employing a list of features that incorporates expression, structure\nand evolutionary conservation. Using the model, they were able to characterize\nevents into “functional” and “non-functional” categories and validate some of\nthe “functional” ones. While there is some excitement regarding the modeling,\nthe curation and assessment of splicing was unclear. However, the\nexperimental follow-up in its current form is challenging to follow as it lacks\ndetails. The paper follows primarily one splicing event in the Tbc1d23-E15 in\nmouse\n1. The paper is dense and written as a short form article, but this makes it very\nchallenging to understand and would benefit from more main figures and a\nsignificant improvement in writing (primarily adding more information).\nResponse: We thank the reviewer for the comments and suggestions. We\nrevised our results section and updated the figures to improve the clarity\n(Revised Fig. 1, Revised Extended Data Fig. 2, 3, 4, 11, 12), per your\nsuggestion.\n2. Experimental sections lack of sufficient information regarding species, model\nsystems used, validation etc. Much of this information is crucial for an\nappropriate review.\nResponse: We apologize for the lack of sufficient information and related\nsections were revised as follows:\n(1) In the ‘Results’ section:\nFor mouse model (Pg. 11, ln. 314-318):\n“To investigate the role of TBC1D23-E15 in vertebrate hematopoiesis, we\ngenerated a transgenic mouse line in which Tbc1d23-E15 was deleted (Fig.\n4a). In contrast to wildtype (WT) mice expressing both long and short\nisoforms of Tbc1d23, the Tbc1d23ΔE15/ΔE15 (E15-KO) mice express only the\nshort isoform (Fig 4b, Extended Data Fig. 7d) with the overall Tbc1d23\nprotein level was comparable to WT mice (Fig. 4c, d)”\nFor zebrafish (Pg.12, ln. 321-326):\n“In addition, to assess the conservation of this phenotype, we targeted the\nlong isoform of Tbc1d23 in zebrafish (Tbc1d23-L, including E15) by\nmorpholino (MO) treatment. The Tbc1d23-L-MO (E15-MO) treatment\nspecifically reduced the mRNA level of the long isoform while increasing the\nexpression of the short isoform, with no significant effect on total Tbc1d23\nprotein level comparable to the fish treated with the control MO (Ctrl-MO)\n(Fig. 4g-i, Extended Data Fig. 7g).”\nFor cell lines (Pg.13, ln.357-359):\n“We generated three in both K562 and HEK293T, one lacking TBC1D23\nexpression (KO) and the other two expressing only TBC1D23-L or\nTBC1D23-S (Extended Data Fig. 9a, Fig. 5c, d)”\n(2) The ‘Methods’ section (Pg.31-32, ln.895-933):\n“Mice. C57BL/6J mice were purchased from GemPharmatech Co. Ltd. The\nTbc1d23ΔE15/ΔE15 (E15-KO) mice were generated in the C57BL/6J\nbackground by GemPharmatech Co. Ltd, whose introns flanking exon 15 of\nTbc1d23 were targeted by 4 sgRNAs using CRISPR/Cas9 technology. E15-\nKO and wildtype littermate controls (WT) of the same sex aged between 8\nand 12 weeks were used in all comparison experiments. Rosa-Cas9 knock-\nin mice were originally purchased from the Jackson Laboratory and\nbackcrossed onto C57BL/6J background. Both heterozygous (Cas9+/-) and\nhomozygous (Cas9+/+) Rosa-Cas9 knock-in mice were used in this study and\ncollectively designated as Cas9+ (Extended Data Fig. 6c). The oligo\nsequences of sgRNA and genotyping primers were listed in Supplementary\nTable 10.\nAll mice were housed and bred under SPF (Specific Pathogen Free)\nconditions at the Laboratory Animal Center of West China Second University\nHospital and were allowed access to diet and water ad libitum. All animal\nexperiments were carried out following the protocols approved by the ethics\ncommittee of West China Second University Hospital [(2018) Animal Ethics\nApproval No.004].\nGenotyping. The Mouse Tail Direct PCR kit (FORE GENE, #TP-01331) was\nused to extract DNA from mouse ails. Genomic DNA was extracted from\nsorted GFP+ cells using the TIANamp Genomic DNA Kit (TIANGEN, #DP304)\nand cells were lysed in elution buffer containing proteinase K and incubated\nat 65 °C for 1 hour, followed by heat inactivation at 95 °C for 5 minutes. PCR\nwas performed for genotyping with primer sequences were listed in\nSupplementary Table 10.\nZebrafish. All zebrafish (Danio rerio) experiments were performed according\nto standard procedures as previously described56,99,100. Both adult fish and\nembryos were maintained at 28.5°C in Aquatic Ecosystems. AB (wild-type)\nand Tg (gata1a:DsRed) strains were used in this study. Larvae were imaged\nat 48 hours postfertilization (hpf), a developmental stage at which sex is not\nyet distinguishable. All experimental protocols were approved by the Animal\nEthical Committee, West China Hospital of Sichuan University.\nMorpholino injection and validation. Antisense morpholino oligonucleotides\n(MO) were purchased from Gene Tools. Tbc1d23-E15-MO (E15-MO) targets\nthe exon 15 of the zebrafish Tbc1d23 gene. Control-MO (Ctrl-MO) is a\nstandard, mismatched control. MO were injected (3 ng per embryo) into the\nembryo at the one-cell stage. To determine the efficiency of MO, total RNAs\nof approximately 40 zebrafish embryos at 48 hpf were extracted and were\nreverse transcribed to complementary DNA (cDNA) for RT-qPCR (qPCR)\nand RT-PCR (see below). Oligo sequences of MO and primers were listed in\nSupplementary Table 10.”\n3. For validation of mechanism, omic-scale experiments are offered such as\nscRNA-seq, IP-mass spec, mass spec for sumo but these are provided with\nlittle context or interpretation.\nResponse: We apologize for not providing sufficient details and updated the\n‘Methods’ and ‘Results’ sections in the revised MS to include relevant\ninformation.\n(1) We added the details of IP-mass spectrometry and SUMO proteome\nanalysis in ‘Methods’ (P g.40-42, ln.1167-1217; Pg.42-43, ln.1235-1270):\n“Immunoprecipitation and Immunoblotting. Cells were harvested and\nlysed in RIPA (Thermo Fischer Scientific, #89900) containing Mini EDTA-free\nprotease inhibitor cocktail (Thermo Fischer Scientific, #A32955) on ice for 30\nmin, and pelleted. The supernatant was solubilized in 0.1 M Tris-HCl\ncontaining 2% SDS and 0.05 M DTT at 95% for 8 min.\nFor immunoprecipitation (IP), cells were harvested and lysed in Western-IP\nlysis buffer (Beyotime, #P0013) containing 10mM protease inhibitors\n(Bimake, #B14002), 20mM N-Ethylmaleimide (Sigma, #128-53-0), 1mM DTT\nsimilar to previous studies(Wang, Cao et al. 2019). Lysates were centrifuged\nat 14,000g at 4 °C for 20 min. Rabbit anti-TBC1D23 (Proteintech, #17002-1-\nAP), mouse anti-FLAG (Proteintech, #20543-1-AP), Rabbit anti-HDCA1\n(CST, #34589) antibodies were equilibrated with Tris-buffered saline (TBS)\n(Sangon Biotech, #C520002) before use. The supernatant from cell lysates\nwas mixed with the protein A/G magnetic beads (Bimake, #B23202) and\nincubated with antibodies at 4 °C overnight. The protein-antibody-magnetic\nbeads complexes were then washed using TBS three times and eluted with\nsodium dodecyl sulfate (SDS) loading buffer (YEASEN, #20315ES) at 95 °C\nfor 10 min. The solutes were loaded and separated by 10% SDS–\npolyacrylamide gel electrophoresis (SDS-PAGE) gel (Epizyme, #PG112) for\nmass spectrometry or immunoblotting.\nProtein extracts or eluted proteins from the IP assay were separated on 10%\nSDS polyacrylamide gels and transferred to 0.45 μm polyvinylidene\ndifluoride (PVDF) membranes (Beyotime, #FFP28). The membranes were\nblocked in TBS buffer with 0.1% Tween 20 and 5% skimmed milk at RT for\n1 hr. The blots were probed with Rabbit anti-TBC1D23 (1:1000), Mouse anti-\nHDCA1 (1:1000) (CST, #5356), Rabbit anti-sumo1 (1:1000) (HUABIO,\n#ET1606-53), Rabbit anti-Sumo2/3 (1:1000) (HUABIO, #ET1701-17),\nMouse anti-Flag (1:1000), Mouse anti-GFP (1:1000), Mouse anti-HA (1:1000)\nor Mouse anti-RanBP2 (1:1000) (Santa Cruz Biotechnology, #sc-74518) at\n4 °C overnight. Then the membranes were incubated with the Goat anti-\nRabit IgG Seconadary antibody HRP conjugated IgG (1:2000) (SAB,\n#L3012-2) or Goat anti-Mouse IgG Seconadary antibody HRP conjugated\nIgG(H+L) (1:2000) (SAB, #L3032-2) at RT for 1 h. The chemiluminescence\nsignal was detected using a ChemiDOCTMMP Imaging System (Bio-Rad\nLaboratories). Image J software was used to analyze the gray value of the\nWB.\nIP-Mass. Gel pieces were destained in 50 mM NH₄HCO₃/50% ACN,\ndehydrated with ACN, and subjected to reduction with 10 mM DTT (56°C, 60\nmin) and alkylation with 55 mM iodoacetamide (room temperature, dark, 45\nmin). After washing, gel pieces were digested with 10 ng/μL trypsin at 37°C\novernight. Peptides were extracted with 50% ACN/5% FA followed by 100%\nACN, dried, and reconstituted in 2% ACN/0.1% FA. The tryptic peptides were\ndissolved in 0.1% formic acid (solvent A), directly loaded onto a home-made\nreversed-phase analytical column (15-cm length, 75 μm i.d.). The gradient\nwas comprised of an increase from 6% to 23% solvent B (0.1% formic acid\nin 98% acetonitrile) over 16 min, 23% to 35% in 8 min and climbing to 80%\nin 3 min then holding at 80% for the last 3 min, all at a constant flow rate of\n400 nl/min on an EASY-nLC 1000 UPLC system. MS analysis was performed\non a Q Exactive™ Plus mass spectrometer with NSI ionization (2.0 kV). Full\nMS scans (m/z 350–1800) were acquired at 70,000 resolutions, followed by\ndata-dependent MS/MS scans (NCE 28, 17,500 resolution) with a dynamic\nexclusion of 15.0 s. Data were processed using Proteome Discoverer search\nengine (v.2.4). Tandem mass spectra were searched against Homo Sapiens\n(SwissProt, 20366 entries) database with trypsin/P digestion (max 2 missed\ncleavages), precursor mass tolerance of 10 ppm, fragment tolerance of 0.02\nDa, carbamidomethyl (C, fixed), and oxidation (M, variable). Peptides were\nfiltered at high confidence with ion score > 20.\n…\nSUMO proteome analysis. TBC1D23 KO K562 cells (4×107), or cells\nexpressing only TBC1D23-L or -S were collected for quantification of SUMO\nmodifications using a PTMScan HS Ubiquitin/SUMO remnant kit (Cell\nSignaling Technologies, #59322) according to the manufacturer’s\ninstructions. In brief, the cells were harvested in 9M Urea Lysis buffer. The\nlysate was sonicated using a sonicator (BioSafer, #650-92) with a microtip at\n20,000 g for 15 min at room temperature and transfer the protein extract\n(supernatant) into a new tube. Protein concentrations were determined by\nbicinchoninic acid (BCA) assay (Thermo Fisher Scientific, #23227). DTT\n(1.25 M) was used for protein reduction, and then iodoacetamide solution\n(Thermo Fisher Scientific, #A39271) was added to the cleared cell\nsupernatant, and diluted 4-fold with 20 mM HEPES pH 8.0 to a final\nconcentration of approximately 2 M urea, 20 mM HEPES, pH 8.0. The\nsamples were digested overnight at 37 °C with WALP enzyme (Cell Signaling\nTechnology, #33036) 0.4 mg/mL stock at the ratio of 1 mg enzyme: 100 mg\nsubstrate. Complete digestion was confirmed by SDS-PAGE. The digested\nsample was acidified with 1/20 volume of 20% trifluoroacetic acid to a final\nconcentration of 1%. The pH was checked by spotting a small amount of\npeptide sample on a pH strip (the pH should be under 3). After acidification,\nallow the precipitate to form by letting the sample stand for 15 min on ice.\nCentrifuge the acidified peptide solution for 15 min at 1,780 g at room\ntemperature to remove any precipitate. The peptide-containing supernatant\nwas transferred into a new 50 mL conical tube without dislodging the\nprecipitated material. Acidified and cleared digest was loaded onto the C18\ncolumn, which was pre-wetted with 5 mL 100% ACN and washed\nsequentially with 1, 3, and 6 mL of Peptide Purification Equilibration Solution\n(0.1% TFA). The column was then washed sequentially with 1, 3, and 6 mL\nof Peptide Purification Equilibration Solution (0.1% TFA) and with 2 mL of\nPeptide Purification Wash buffer. Place columns above new 50 mL\npolypropylene tubes to collect eluate. Elute peptides with a sequential wash\nof 3 mL and then 7 mL of Peptide Purification Elution Solution (0.1% TFA,\n50% acetonitrile). TFA was removed by freezing the eluate on dry ice (or -\n80°C freezer) overnight and lyophilizing frozen peptide solution for a\nminimum of 2 days. The tube containing lyophilized peptides was centrifuged\nfor 5 minutes at 2,000 g at room temperature to collect material for dilution\nin HS IAP Bind buffer #1. PTMScan HS antibody-bead slurry in ice-cold\nphosphate-buffered saline (PBS) was added to each clarified sample.\nSamples were incubated on an end-over-end rotator at 4 °C for 2 h. The\nbeads were collected using a magnetic stand and then washed sequentially\nwith HS IAP Wash buffer and LCMS water to wash the beads. The beads\nwere resuspended in IAP Elution buffer (0.15% TFA) and the eluted sample\nwas collected in a new microcentrifuge tube.”\nAdditionally, we added 2 schematic diagrams to illustrate the experimental\ndesigns and workflow for both IP-mass spectrometry and SUMO proteome\nanalysis (Revised Fig. 5a, g).\nRevised Fig. 5a, Scheme of IP-MS study design created with Biorender.com.\nPlasmids encoding empty mCherry vector (mCherry), mCherry-TBCD23-L, or\nmCherry-TBCD23-S were overexpressed in TBC1D23-KO K562 cells, followed\nby IP assays using an anti-TBC1D23 antibody.\nRevised Fig. 5g, Diagram of the experimental design and workflow for MS analysis of\nSUMOylation peptides affected by TBC1D23-L or TBC1D23-S. The plasmids of\nmCherry, mCherry-TBCD23-L and mCherry-TBC1D23-S were overexpressed in\nTBC1D23-KO K562 cells. Cells were lysed in a urea-based buffer, followed by WALP\ndigestion of the proteins. The peptides were purified by reversed-phase, solid-phase\nextraction. Using magnetic antibody beads targeting the SUMO remnant motif (K-ε-\nGG), SUMOylation peptides were then captured by immunoprecipitation and eluted to\nconcentrate them for LC-MS/MS analysis. Created with BioRender.com.\n4. While I appreciate the work for obtaining FAScores, I’m not sure that the\nexperimental components and the selection of AS for validation necessarily\nrequires this score. One might prioritize the same events based on conservation,\nprotein domain impact, and lineage. Thus, more should be done to discuss why\nthe model facilitates discovery beyond the standard practices of looking at\ndynamics and conservation.\nResponse: Thanks for the reviewer’s comments. We would like to clarify how\nFAScore aides in the identification of functional AS in this study.\nFirstly, while sequence conservation and protein domain are among the most\nimportant features, they are not the only determinants for prediction. Our model\nintegrates 19 features into a unified FAScore and classified the candidates into\n‘functional AS’, ‘non-functional AS’, and ‘uncertain AS’. We observed the\ncorrelation between conservation and FAScore was not as high as one may\nimagine (Revised Extended Data Fig. 3a). Notably, FAScore successfully\nidentifies many functional AS events with low sequence conservation. For\ninstance, known functional AS events such as Kdm1a-E7 (Zibetti, Adamo et al.\n2010) and EHMT2-E10 (Fiszbein, Giono et al. 2016) exhibit low conservation\nbut high FAScores (Revised Extended Data Fig. 3b). This demonstrates that\nrelying solely on conservation for prioritization would miss these important\ncandidates.\nSecondly, the features used in the model were of fine granularity, for example,\nfor the protein domain features, we interrogated features including structural\nhomologs and integrity, domain integrity, transmembrane helices, signal\npeptides, and target peptides. Our model can quantify the contribution of\ndifferent features (Revised Extended Data Fig. 2d), avoiding the subjective\nbiases in rule-based filtering approaches.\nFinally, unlike traditional filtering methods that yield a large set of candidate AS\nevents without internal prioritization, FAScore provides a continuous score. This\nenables ranking of candidates by their relative score, thereby offering a\nquantitative framework for prioritizing targets and increasing the efficiency of\ndownstream experimental validation.\nWe also added text in the discussion regarding why the model facilitates\ndiscovery beyond the standard practices of looking at dynamics and\nconservation. (Pg.16, ln. 448-454)\n“Our model integrates multiple features into a unified predictive score, and\ncan automatically quantify the contribution of different features to the\nprediction, avoiding the subjective biases in features-based filtering\napproaches. Finally, compared to traditional selecting methods that yield a\nlarge set of candidate AS events without prioritization, FAScore provides a\ncontinuous prediction score. This enables ranking candidates by their\nrelative score, thereby offering a quantitative framework for prioritizing\ntargets and increasing the efficiency of downstream functional validation.”\nRevised Extended Data Fig. 3a, Density plot showing the correlation between\nFAScore and evolutionary conservation across species in ES events, with darker\nregions indicating higher density of events. The red circle highlights a region containing\nsplicing events with high FAScore but low conservation. Spearman correlation\ncoefficient and corresponding p-value are shown.\nRevised Extended Data Fig. 3b, Scatter plot showing the correlation between\nFAScore and evolutionary conservation score across species. Each point represents\na known functional splicing event.\nRevised Extended Data Fig. 2d. Importance ranking of features used in the final\nmodel by average Gini index\n5. Description of “Dynamically expressed genes” is difficulty to follow. This\nsection needs rewriting and perhaps schematics to help explain the method in\nmore detail. Representative example data for certain genes or alternative\nsplicing events would be helpful.\nResponse: We thank the reviewer for this constructive suggestion. We\nincluded a schematic diagram to visually illustrate the workflow for identification\nof dynamically expressed genes and AS events across hematopoietic lineages\nor developmental stages (Revised Fig. 1b). We also incorporated examples of\ndynamically-expressed gene and splicing events, such as MEIS1, MEIS1-E12,\na known functional AS event (Zeddies, Jansen et al. 2014), and MEIS1-E1\n(Revised Fig. 1c-e).\n(Pg. 5-6, ln. 124-160)\n“Genes or AS events that are dynamically regulated during developmental\nor cell differentiation trajectories are more likely to play critical roles in these\nbiological processes, as their temporal regulation often reflects direct\nfunctional involvement in key cellular transitions or fate decisions34. To\nquantify the dynamic changes in gene expression and alternative splicing,\nwe computed a series of feature scores based on the expression levels or\nPSI values across developmental time points or differentiated cell types (Fig.\n1b). Specifically, we calculated the expression breath (Gene.Range and\nAS.Range) to assess the magnitude of change, the correlation coefficient\n(Gene.cor and AS.cor) and its statistical significance (Gene.cor.p and\nAS.cor.p) to evaluate the association with temporal or differentiation\ntrajectories, the linear regression slope (Gene.slope and AS.slope) and its\nstatistical significance (Gene.slope.p and AS.slope.p) to determine the rate\nand statistical significance of the progressive changes, as well as tau\n(Gene.tau and AS.tau) to measure expression specificity (Fig. 1b, see\nMethods). As a representative example, the MEIS1, a transcription factor\nwhose inactivated impaired hematopoietic stem cell niche development,\nmegakaryocyte and platelet deficiencies in mice42-44. We quantified six\nfeatures of the dynamic gene expression of MEIS1. It showed a gradually\ndecreased expression in hematopoiesis with an upregulation in\nmegakaryocyte (Fig. 1c), consistent with a previous report45. Meanwhile,\nMEIS1 contains two ES events, the exon (E) 1 and 12. The MEIS1-E12\nexhibits high level of changes (Fig. 1d), consistent with its function in\ninducing the fate decision of human hematopoietic progenitors toward\nmegakaryocyte-erythroid progenitors (MEPs)46. In contrast, the MEIS1-E1\nshowed a stable expression across cell types with low values of dynamic\nfeatures (Fig. 1e).\nTo facilitate intuitive assessment of the dynamics of genes and AS events,\nwe integrated the dynamic features described above into a unified dynamism\nscore (DyScore). Based on its bimodal distribution, we applied Gaussian\nMixture Modeling (GMM) to classified the genes or AS events into \"dynamic\"\nand \"non-dynamic\" (Extended Data Fig. 1c, e, see Methods). Overall, 27%\nof AS events and 56% of genes were classified as dynamic (Extended Data\nFig. 1d, f). Notably, in both human and mouse, dynamically regulated AS\nevents were more likely to be cell lineage-specific, with 63% dynamic in only\none hematopoietic lineage compared to 29% of dynamically expressed\ngenes (Extended Data Fig. 1g). Moreover, dynamically regulated AS\nexhibited significantly higher evolutionary conservation scores in the regions\nflanking the splice sites across all lineages and fetal hematopoietic organ\n(paired t-test, p <1.5x10-4, Extended Data Fig. 1h). These results suggest\nthat dynamically regulated AS events may play specialized and conserved\nregulatory roles in hematopoietic lineage, highlighting the importance of\nsplicing dynamics in the spatiotemporal control of transcriptome complexity\nduring hematopoiesis.”\nRevised Fig.1. The examples of dynamically regulated gene and AS events.\nb, Schematic representation of dynamic profiling for gene expression and alternative\nsplicing (AS). A linear fit (blue dashed line) and a Loess curve (green solid line) are\nshown, with the shaded area representing the 95% confidence interval. Derived\nfeatures include: range of genes (Gene.range) or AS (AS.range) to assess variation\nbreadth; correlation coefficient (Gene.cor, AS.cor) and its p-value (Gene.cor.p,\nAS.cor.p) to quantify association with temporal or differentiation trajectories; linear\nregression slope (Gene.slope, AS.slope) and its p-value (Gene.slope.p, AS.slope.p) to\nestimate rate and significance of change; and tau (Gene.tau, AS.tau) to evaluate\nspecificity.\nc-e, Expression and splicing patterns of MEIS1 and its two AS events during erythroid\ndifferentiation. (c) Expression values (transformed TPM) of MEIS1 (Left) and the\nvalues of dynamic features (Right). Percent spliced in (PSI) values (Left) and dynamic\nfeature scores (Right) for the exon 12 skipping event in transcript MEIS1-201 (d) and\nthe exon 1 skipping event in transcript MEIS1-202 (e). A linear fit (blue dashed line) and a\nLoess curve (green solid line) are shown, with the shaded area representing the\n95% confidence interval.\n6. For expression breadth, it is mentioned that TPMs for each gene was used.\nTPMs are normalized within samples and are not the best measure for\nbetween-sample comparisons. Scaling it between 0 and 1 could lead to\nconfounding effects. What does the overall distribution of dynamic scores look\nlike for genes and AS events. How was 0.4 chosen?\nResponse: We used TPMs because TPMs normalize both sequencing depth\nand gene length(Zhao, Li et al. 2021). This double normalization makes TPM\nvalues comparable both between samples (accounting for sequencing depth)\nand between genes (accounting for gene length) within a sample. This\ncomparability is essential for robust calculation of dynamic scores and for\nmeaningful cross-gene comparisons (such as expression range) of the scores.\nMethods such as DESeq2 (based on library size) (Love, Huber et al. 2014) are\ndesigned for comparisons between samples/conditions; however, it is not ideal\nfor comparisons across genes needed for dynamic scores since DESeq2 does\nnot normalize gene length.\nIn a comparison analysis, we found high correlations between TPM and\nDESeq2 normalized counts within the same samples (median Spearman ρ =\n0.96), and we also found high consistency between two normalisation methods\nfor the same gene (median Spearman ρ = 0.80), such as Paqr8, Eya1 and\nStau2. Therefore, we maintained the use of TPM as it inherently provides the\ncross-gene and cross-sample comparability which is fundamental to our\nanalysis of the expression breadth and dynamics. The additional analysis was\nin included in Revised Extended Data Fig. 11a-c.\nSecond, the expression breadth for each gene (Gene.range) was calculated by\ncomputing the difference (max - min) of the log2-transformed values (log2(TPM\n+ 1)) within each gene (Revised Extended Data Fig. 11d). Therefore, when the\nexpression breath is 0, it means no difference. Meanwhile, when the expression\nbreath is 1 or greater than 1, it means a minimum 2-fold expression difference\n(log2(FC) = 1) within a gene. We use this cutoff to category a gene as highly\nvariable (Revised Extended Data Fig. 11e). The rationale is that the DyScore\nof a gene is calculated by adding these features together, so we need the\nexpression breadth to have an upper limit. Otherwise, the DyScore may be\ninflated by a single feature. For AS events, we directly computed the difference\nof PSI values within each AS event (AS.range). Since PSI is intrinsically\nbounded between 0 and 1, the AS.range values also naturally fall within the 0\nto 1 range, consistent with genes (Revised Extended Data Fig. 11f).\nThird, we added the distribution of absolute dynamic score (Revised Extended\nData Fig. 1c, e). Based on its bimodal distribution, we applied Gaussian\nMixture Modeling (GMM) to classify the genes or AS events into \"dynamic\"\nand \"non-dynamic\". A GMM with k=2 components effectively captured the\nbimodal distribution (G1 and G2). The optimal threshold separating the two\ncomponents was determined by the probability density intersection method,\nresulting in a cut-off of 0.467 for genes (Revised Extended Data Fig. 1e) and\n0.324 for AS events (Revised Extended Data Fig. 1c). We then calculated the\narithmetic mean of the gene-level and AS-level cutoffs, which was 0.40.\nWe revised this section in ‘Methods’ as follows:\n(1) Expression breadth: range value (R) (Pg.19, ln.548-559)\n“The value range (R) represents the difference between the maximum and\nminimum expression values, reflecting the extent of variation. For genes,\nexpression breadth (R ) was calculated as follows: a pseudo-count of 1\ngene\nwas added to the TPM values, which were then log2-transformed; R was\ngene\nthen computed as the difference (max − min) among the transformed values\nfor each gene (Extended Data Fig. 11d). when the expression breath is 0,\nit means no difference. Meanwhile, when expression breath is 1 or greater\nthan 1, it means a minimum 2-fold expression difference (log2(FC) = 1) within\na gene. Thus, we applied a transformation step wherein any R value ≥ 1\ngene\nwas set to 1 (Extended Data Fig. 11e). For AS events quantified using\nPercent Spliced In (PSI), the range (R ) was directly computed as the\nAS\ndifference in PSI values for each AS event. Since PSI is naturally bounded\nbetween 0 and 1, RAS values also inherently fall within the 0–1 range\n(Extended Data Fig. 11f).”\n(2) Definition of dynamic genes and AS (Pg.22, ln.620-628)\n“Based on bimodal distribution of absolute dynamic score, we applied\nGaussian Mixture Modeling (GMM) to classify the genes or AS events into\n\"dynamic\" and \"non-dynamic\". A GMM with k = 2 components effectively\ncaptured the bimodal distribution. The optimal threshold separating the two\ncomponents was determined by probability density intersection method,\nresulting in a cut-off of 0.467 for genes and 0.324 for AS events. We then\ncalculated the arithmetic mean of the gene-level and AS-level cutoffs, and it\nwas 0.40. Finally, genes or AS events were divided into two classes,\nincluding dynamically expressed genes/AS events (>= 0.4 or <= -0.4), and\nnon-dynamically expressed genes/AS events (-0.4~0.4).”\nRevised Extended Data Fig. 11. The distribution of expression breadth of genes\nand alternative splicing events.\na-c, Scatter plots showing the correlation between TPM and DESeq2-normalized\nexpression values for three example genes: (a) Paqr8 (R = 0.94, p < 2.2×10-16), (b)\nEya1 (R = 0.96, p < 2.2×10-16), and (c) Stau2 (R = 0.96, p < 2.2×10-16). The blue line\nin each panel represents the linear regression fit.\nd-f, Density distributions of expression breadth of genes and AS, including Gene.range\n(d), Capped Gene.range (e) and AS.range (f).\nRevised Extended Data Fig. 1. c, e, Distribution of absolute dynamic score (DyScore)\nshown as overlaid histogram and KDE for AS events (c) and genes (e). G1 and G2\ndenote the two fitted Gaussian components used for classifying dynamic and non-\ndynamic elements.\n7. For expression gradient, it’s important to show that the relationship was linear\nbetween time points and gene expression and AS events which might not\nalways be the case. How were these dealt with? Furthermore, no reasoning\nwas provided for capping a continuous scale such as regression coefficient to -\n1 and 1. Examples of events that agree with expression specificity score would\nbe helpful.\nResponse: We agree that temporal relationships between time points and\ngene expression and AS events are not necessarily linear. We used generalised\nadditive modeling (GAM) and effective degrees of freedom (edf) to objectively\ncharacterize relationship linearity. In GAMs, the edf is a measure of a smooth\nterm's \"wiggliness\" or complexity, reflecting how much the model is allowed to\ndeviate from a straight line. According to the common classification method of\nedf(Hunsicker, Kappel et al. 2016), AS or genes can be divided into three\ncategories: linear (edf ≈ 1), weakly non-linear (1 < edf < 2) and strongly non-\nlinear (edf > 2) (Revised Extended Data Fig. 12a). The result showed\npredominant linear tendencies (51.8%), with edf values clustering near 1\n(Revised Extended Data Fig. 12a).\nAdditionally, we found slope value from linear regression was representative in\nthe three categories. For linear patterns (edf ≈ 1), slopes provided optimal\nquantification of the change rates (e.g., MED24; Revised Extended Data Fig.\n12b). For weakly nonlinear patterns (1 < edf < 2), slopes remained suitable\nchange rate indicators in 51% of cases (validated by >50% variance explained)\n(Revised Extended Data Fig. 12b, e.g., MAD1L1). Within strongly nonlinear\n(edf > 2), slopes effectively captured the change rates for genes/AS exhibiting\nmonotonic trends (Revised Extended Data Fig. 12b, e.g., FGG). However, for\nstrongly nonlinear features with non-monotonic oscillations, the sign of slopes\nrepresents only the direction of change (Revised Extended Data Fig. 12c, e.g.,\nHCCS). To address this, we introduced Range metrics as the complementary\nglobal descriptors. Range provides supplementation when slopes\nunderestimate dynamics. Thus, linear fitting is a important feature we assess\nfor dynamic gene or AS regulation. In sum, the majority of temporal\nrelationships between time and expression and AS events are linear or close to\nlinear. To simplify the process and comparability, we used linear regression in\nour study.\nFor the regression coefficient (slope), to address potential skewness in gene\nexpression distributions while preserving biological fidelity, we performed the\nlog (TPM+1) transformation for gene expression values followed by linear\n2\nregression against normalized timepoints. While this transformation mitigated\nbiases from highly expressed genes, initial slopes spanned from -2 to 2 (with\nthe majority of them distributed around 0). However, these slopes exhibited only\nlow correlation with raw TPM slopes (Pearson r= 0.17, Revised Extended Data\nFig. 12d). Subsequently, we applied the capping at [-1,1] to extreme values\n(affecting <1% of genes), and this optimization substantially improved the\ncorrelation to r of 0.84 (Revised Extended Data Fig. 12e), suggesting an\nimproved biological relevance. The optimized slope values also enabled direct\nadditive integration of different features, providing a mathematically consistent\nfoundation for downstream dynamic score calculation. For AS events, PSI\nvalues inherently generated slopes within [-1,1] without processing (Revised\nExtended Data Fig. 12f), maintaining analytical consistency with the\nprocessed gene expression metrics.\nWe revised the ‘Methods’ section (Pg.20-21, ln.576-600) as follows:\n“To objectively characterize the linearity of relationships, we employed\ngeneralized additive modeling (GAM) and effective degrees of freedom (edf),\nwhich quantifies the complexity or “wiggliness” of a smooth term, reflecting\nthe extent to which the model deviates from a straight line. The observation\nthat edf values clustered near 1(Hunsicker, Kappel et al. 2016) indicates that\ngene expression and splicing levels exhibit predominantly linear tendencies\n(51.8%) (Extended Data Fig. 12a-c). Accordingly, we primarily conducted\nlinear regression analyses, where the regression coefficient (β) represents\nthe magnitude of influence exerted by the independent variable on the\ndependent variable, and the p-value of the regression variable (P ) reflects\nβ\nthe statistical significance of this association.\nTo mitigate potential skewness in gene expression distributions while\npreserving biological interpretability, we applied a log2(TPM + 1)\ntransformation to gene expression values, which were then subjected to\nlinear regression against normalized differentiation cell types or\ndevelopmental time points. Although this transformation reduced bias from\nhighly expressed genes, the resulting slopes initially ranged from -2 to 2, with\nmost values concentrated near 0. These initial slopes showed only weak\ncorrelation with slopes derived from raw TPM values (Pearson r = 0.17,\nExtended Data Fig. 12d). We therefore applied capping to extreme values\nat the range [-1, 1], which affected fewer than 1% of genes. This adjustment\nsubstantially improved the correlation to r = 0.84 (Extended Data Fig. 12e),\nsuggesting enhanced biological relevance. The optimized slope values also\nsupport direct additive integration across different features, providing a\nmathematically consistent basis for downstream dynamic score calculation.\nFor alternative splicing (AS) events, PSI values inherently yield slopes within\n[-1, 1] without additional processing (Extended Data Fig. 12f), ensuring\nanalytical consistency with processed gene expression metrics. Furthermore,\np-values were binarized: values below 0.05 were assigned a significance\nlabel of 1, and all others were labeled 0.”\nRevised Extended Data Fig. 12. Systematic validation of linear modeling\napproach as dynamics features.\na, Distribution of effective degrees of freedom (edf) from generalized additive models\n(GAM). Categories: linear (edf ≈ 1.0), weakly nonlinear (1 < edf < 2), strongly\nnonlinear (edf > 2).\nb, Dynamics of expression of representative genes (MED24, MAD1L1, FGG) across\nedf pattern categories. The yellow solid lines represent the GAM fits, while the blue\ndashed line indicates the linear regression fits. Gray shaded areas represent 95%\nconfidence intervals.\nc, Dynamics of expression of representative genes (HCCS) in high edf and near-zero\nslope. The yellow solid lines represent the GAM fits, while the blue dashed line\nindicates the linear regression fits. Gray shaded areas represent 95%\nconfidence intervals.\nd, Correlation between expression slopes of raw and log-transformed TPM expression\nper gene. Pearson's correlation coefficients (R) and p-value are labeled.\ne, Correlation of capped slope between raw and log-transformed TPM expression per\ngene. Pearson's correlation coefficients (R) and p-value are labeled.\nf, Density distribution of slopes of alternative splicing based on PSI.\n8. Alternative splicing patterns are already shown to be lineage specific (Merkin\nof that.\nResponse: Merkin et al. (2013 Science) used \"lineage\" in the context of\nevolutionary divergence (e.g., rodent-specific, primate-specific, or mammalian-\nconserved AS events). This analysis compared AS patterns across species\n(mouse, rat, macaque, cow, chicken) to identify phylogenetically restricted\nsplicing and established the concept of lineage-specific alternative splicing (AS)\nin evolutionary biology. The word \"lineage\" in our study refers to the\nhematopoietic cell lineage during including myeloid, erythroid and lymphoid\nlineages. To avoid confusion, we used the wording of “cell lineage-specific” in\nthe revised MS.\n9. Sequence conservation is part of the formula used to calculate the score and\nis one of the top features based on gini index. Could the model be biased\ntowards events with high conservation (i.e. is this circular)? Conservation and\nprotein domain composition is a key feature of gene importance, thus is an ML\napproach required for identifying these events and prioritizing them?\nResponse: Although sequence conservation and protein domain are the top\ncontributing features, they are not always the decisive ones, given that our\nmodel integrates multiple features into a unified predictive score.\nAs shown in the Revised Extended Data Fig. 3a, the conservation and\nFAScores are weakly correlated, and FAScores predicts many AS events with\nlow sequence conservation. Indeed, when we plotted the conservation and\nFAScores of the AS events with known function (published in literature), some\nare of low sequence conservation, such as Kdm1a-E7 (Zibetti, Adamo et al.\n2010), EHMT2-E10 (Fiszbein, Giono et al. 2016) (Revised Extended Data Fig.\n3b). Therefore, if we only use conservation to prioritize AS events, these events\nwould be missed.\nTo address the concerns on the model biased towards events with high\nconservation, the two heat maps in revised Extended Data Fig. 3c, d\ndemonstrated that the conservation and FAScore are weakly related as there\nare events with high conservation but having low FAScore and events with high\nFAScore with low conservation (boxed in revised Extended Data Fig. 3c,d).\nAs to protein domain composition, it ranked the fifth contributing 5.6%\nimportance using Gini index. Moreover, the percentage of functional AS events\nwith defined domain was only 28.5%. The selection based on domain\ncomposition may also miss interesting AS events. For example, in this study we\nshowed that TBC1D23-E15, which only encodes 15 amino acids with no\nindicated domain or motif, plays important roles in erythropoiesis.\nRevised Extended Data Fig. 3. Correlation between FAScore and evolutionary\nconservation across alternative splicing events.\na, Density plot showing the correlation between FAScore and evolutionary\nconservation across species in ES events, with darker regions indicating higher density\nof events. The red circle highlights a region containing splicing events with high\nFAScore but low conservation. Spearman correlation coefficient and corresponding p-\nvalue are shown.\nb, Scatter plot showing the correlation between FAScore and evolutionary\nconservation score across species. Each point represents a known functional splicing\nevent.\nc, Heatmap displays normalized feature scores for high conservation with high\nFAScore (50 were randomly selected) or low FAScore (50 were randomly selected) AS\nevents identified via PU-learning. Each row represents an AS event, and each column\ncorresponds to a predictive feature, with the final FAScore.\nd, Heatmap displays normalized feature scores for high FAScore with high\nconservation (50 were randomly selected) or low conservation (50 were randomly\nselected) AS events identified via PU-learning. Each row represents an AS event, and\neach column corresponds to a predictive feature, with the final FAScore.\n10. Does any non-functional or uncertain event have literature on their\nfunctional relevance? What is the false negatives on this model?\nResponse: Our final model was trained using reliable positive and negative\nsets inferred by the positive set from the literature through a PU-learning\nprocess, with 20% of the data held out as a validation set. The trained model\nachieved a low false negative rate of 0.050 on this validation set (10-fold cross\nvalidation, sd = 0.012).\nAs shown in Figure R1.2, the true positive set in the validation set mostly have\nhigh FAScore. As the publications on functional AS were still limited, we did not\nfind ‘false negative’ examples, however, Mbnl1-E5, was classified as “uncertain”\nwith a score of 0.470, while it was reported to play a functional role during\nterminal erythroid differentiation in mice(Cheng, Shi et al. 2014).\nFigure R1.2 Rank-ordered distribution of FAScore for all AS events in the validation\nset from PU-learning strategy. AS events are ranked by descending FAScore (0.00 to\n1.00). The functional AS events from literature are explicitly labeled.\n11. The statement: \"In short, alternative splicing of TBC1D23-E15 leads to\ndifferential binding affinity to the RanBP2/RanGAP1 complex and may affect\ncellular SUMOylation level.\" is not fully supported by experiment. The raw gels\nindicate minimal difference in interaction and exclusively in an over-expression\nsetting, further the term affinity is used but not measured using biophysical\nassays. Similarly, the impact of sumoylation needs to be further validated\nbeyond Co-Ip experiments.\nResponse: We endeavored to obtain recombinant proteins for RanBP2 and\nRanGAP1 for assays such as SPR, but unfortunately not succeeded. We\ntherefore revised the MS to avoid statement of ‘binding affinity’ and rephased\nthe wording as “binding capacity”.\nTo further elucidate the differential binding capacity of the two isoforms, we\nrepeated the co-IP assay using cell lines expressing only TBC1D23-L or -S\nisoform which was generated by reconstitution of the isoform expression in the\nTBC1D23 full gene KO cell line and the total protein level of TBC1D23 (L or S)\nwas comparable to the control cells (Revised Fig. 5c, d, see below).\nUsing\nthese cell lines, we showed that the level of RanBP2 and RanGAP1 bound to\nTBC1D23-S was 4 and 5 times higher than TBC1D23-L, respectively (Revised\nFig. 5e, f). The additional experiment and the accompanying figures were\nincluded in the revised MS (Pg. 13, ln. 356-368)\n“To investigate how TBC1D23-E15 regulates erythropoiesis, we analyzed\nthe interactomes of long and short isoforms of TBC1D23 by an integrated\nmulti-omics approach. We generated three in both K562 and HEK293T, one\nlacking TBC1D23 expression (KO) and the other two expressing only\nTBC1D23-L or TBC1D23-S (Extended Data Fig. 9a, Fig 5c, d), and\nperformed IP-MS in the K562 cell lines (Fig. 5a, Supplementary Table 6).\nRanBP2 was identified as one of the candidate proteins showing different\nbinding toward TBC1D23-S or TBC1D23-L (Extended Data Fig.9b, c).\nRanBP2 is a known Sumo E3 ligase that forms a complex with SUMO1-\nmodified RanGAP161. The immunofluorescence analysis showed that\nTBC1D23 co-localized with RanBP2 (Fig 5b). Using the TBC1D23-KO,\nTBC1D23-L and TBC1D23-S KEK293T cell lines overexpressing RanBP2 or\nRanGAP1 E3 SUMO ligase domain, we confirmed TBC1D23-S consistently\nprecipitated more RanBP2 and RanGAP1 than TBC1D23-L by 5 and 4 times,\nrespectively (Fig 5e, f) suggesting that E15 exclusion enhances engagement\nwith the complex (Fig 5e, f).”\nAdditionally, the SUMO proteomics analysis, presented in revised Fig. 5i,\nidentified 11 hematopoietic regulators with significantly differential SUMOylation\npatterns in TBC1D23-L and -S expressing cells, including HDAC1-K476.\nFinally, we performed additional Co-IP analysis of SUMOylated HDAC1 using\nbone marrow cells isolated from both WT and E15-KO mice and the results\nwere included in the revised Fig. 5m, n and MS (Pg. 14, ln. 401-402):\n“Furthermore, HDAC1 SUMOylation was elevated by 1.5-fold in BMCs from\nE15-KO mice compared to those of WT controls”\nRevised Fig 5c, Western blot analysis of TBC1D23 protein expression levels in\nHEK293T cells under four conditions: sgRNA-control cells, TBC1D23-KO cells\nreconstituted with either empty Flag vector (Flag), Flag-tagged short isoform of\nTBC1D23 (Flag-TBC1D23-S), or Flag-tagged long isoform of TBC1D23 (Flag-\nTBC1D23-L). Representative results from three independent experiments are shown.\nRevised Fig 5d, Statistical analysis of Tbc1d23 expression level. The level was\ndetermined by normalizing the image gray values in (c) (n=3).\nRevised Fig 5e, Co-immunoprecipitation (co-IP) analysis of the binding capacity of\nTBC1D23-S and TBC1D23-L to the RanBP2 complex. HA-RanGAP1 and EGFP-\nRanBP2-E3 were co-transfected into TBC1D23-KO HEK293T cells reconstituted with\nFlag, Flag-TBC1D23-S, or Flag-TBC1D23-L. Cell lysates were subjected to\nimmunoprecipitation (IP) with anti-Flag antibody and immunoblot with anti-Flag, anti-\nHA, and anti-EGFP antibodies.\nRevised Fig 5f, Images of membranes in (e) were analyzed using ImageJ software to\ncalculate the relative intensity of bands for HA-RanGAP1 and EGFP-RanBP2-E3 in IP\nsample (n=3).\nRevised Fig 5i, Heatmap displaying the abundance of proteins involved in the\nregulation of the hematopoiesis pathway. Abundance values are scaled, with row\nlabels indicating protein names and SUMO-modified sites.\nRevised Fig 5m, Co-IP analysis of the HDAC1 SUMOylation level in WT and E15-KO\nmice. Bone marrow cells (1×10⁸) were lysed and subjected to IP with anti-HDAC1\nantibody and immunoblot with anti-SUMO1 and anti-Hdac1 antibodies.\nRevised Fig 5n, Membrane images in (m) were analyzed using ImageJ software to\nquantify the relative HDAC1 SUMOylation level in E15-KO mice compared to WT mice\n(n=3).\n12. Methods are lacking as is model validation for zebrafish and mouse.\nResponse: We apologize for the lack of information. As shown above in point\n2, we updated the ‘Methods’ section in the revised MS to include the validation\nmethods.\nReviewer #2:\nThis manuscript describes a Random Forest based approach to score\nalternative splicing (AS) events based on their likelihood to have functional\nrelevance for specific hematopoietic lineages. Based on this score (FAScore)\nas well as the differential regulation of sequence inclusion, they then selected\nfour AS events, which they functionally validated. They follow up on one of them\n(exon 15 of TBC1D23), and elucidated the mechanistic details on how the\nskipping/inclusion of that exon leads to its physiological impact.Overall, I found\nthe experimental tests and validations to be very well executed, specially for\nTBC1D23. The development of the FAScore is not particularly original, perhaps\nwith the exception of making it lineage-specific (instead of \"global\"). However,\nthe authors seem to then average FAScores and used them in combination with\nthe DyScore (i.e. classic delta PSI) to select the events of interest (Fig. 2d).\n1. The positive training set to develop the FAScore is of 63 events according to\nSup Table 3 (even though the Methods section says 71). This number seems\nquite low for an RF approach. Also, looking at the references in Table S3, for\nsome of these events I am not sure there is actual evidence of functionality,\nonly of protein translation by Western blot. This should be clarified. Finally, most\nof the events used for training are of the exon skipping type (see point 2 below).\nThese limitations should be explicitly acknowledged in the main text.\nResponse: Thanks for the reviewer’s constructive suggestions. We fully agree\nwith comment regarding the training of exon-skipping (ES) type and the events\nwith ‘actual evidence’. We therefore retained only the exon skipping events in\nthe analysis and the positive events with experimental evidence of functionality,\nwhile others with only Western blot evidence were removed (Kalsotra and\nCooper 2011, Barbosa-Morais, Irimia et al. 2012, Merkin, Russell et al. 2012,\nMazin, Khaitovich et al. 2021) (Revised Supplementary Table 3). Although the\nexon skipping event is the most experimentally studied in the field of alternative\nsplicing, we have tried our best to include ES events with published functions\nand the resulting positive training set consists of 58 ES events with function.\nWe acknowledge that the size of our positive set remains limited for\nconventional machine-learning methods. To address this, we have employed a\nPU-learning (Positive-Unlabeled Learning) framework (Hao, Colak et al. 2015,\nBekker and Davis 2020), which is particularly suited for scenarios where\nconfirmed positive samples are scarce and unlabeled samples are abundant\n(lacking definitive negative labels). This approach is well-aligned with the\nconstraints of our dataset. By the method, we defined a reliable positive set\n(4,803) and a negative set (5,245) to form the final training set. The clear\nseparation between positive and negative sets across most features showed\nthe robustness of the PU-learning selection and the discriminative power of the\nfeature set (Figure R2.1). Moreover, the size of the training set was acceptable\nfor the Random Forest model.\nFigure R2.1. The heat map displays normalized feature scores for 4,803 reliable\npositive and 5,245 reliable negative AS events identified via PU-learning. Each row\nrepresents an AS event, and each column corresponds to a predictive feature, with the\nfinal FAScore. Values are z-score normalized.\n2. It is unclear to me how the benchmarking results in Figures 1 and 2 and in\nthe associated Supplementary Figures look for each type of AS event\nseparately. All events seem to be analyzed and treated together, which I believe\nis potentially problematic. As mentioned above, FAScore has been mainly\ntrained for exon skipping events and it is known that the level of conservation,\nfunctionality, etc. is very different for different types of events.\nResponse: We have updated the positive set of exon-skipping events by\nincorporating functional annotations from the literature, applied PU-learning to\nexpand the positive set and infer the negative set, and reconstructed a\nprediction model specifically for exon skipping events. This updated model\ndemonstrates an improved performance (ROC-AUC 0.993) compared to the\nprevious one (0.931) (Revised Fig. 2a).\nRevised Fig. 2a. Evaluation of FAScore prediction performance using two\ncomplementary metrics, including the Receiver Operating Characteristic (ROC) curve.\nThe grey diagonal line represents the performance of a random classifier. AUCmax:\nthe max area under the curve (AUC) value in 10-fold-validation.\n3. Another aspect I found unclear is how the authors have dealt with constitutive,\ncryptic and \"real\" AS exons. Based on the number of AS events analyzed, it\nseems to me they have applied their score to all exons in the genome. What\ndoes it mean if a constitutive exon has a high FAScore? Does it mean that its\nmis-splicing will have a functional impact (as expected)? Again, it would be\ngood to see how the FAScore ranges for the different types of events\n(constitutive, AS, etc.), and how the benchmarkings perform separately per\ninclusion type. For instance, if all highly regulated exons (high DyScore) have\nhigh FAScores, this means the FAScore is mainly a non-useful score on its own.\nThis score seems to categorize thousands of exons/events as functional, which\nis likely questionable.\nResponse: In our study, we focused only on alternative splicing exons without\nconstitutive exons. Our updated methods were trained based on exon-skipping\nevents. We do not think that it can infer the FAScores for constitutive exons,\ngiven the level of conservation and genomic features between constitutive and\nAS exons are different.\nOur study includes annotated (\"real\") and unannotated (“cryptic/novel”) exons,\nto minimize the potential of noise, we included splice junctions that meet the\nfollowing criteria: (1) independent detection by ≥10 uniquely mapped reads in\nat least two individual samples, (2) ≥100 uniquely mapped reads supporting the\njunction when summed across all samples, and (3) unannotated splice\njunctions required to possess one of the canonical intron motifs (GT/AG,\nGC/AG, or AT/AC). Next, we compared the FAScores of annotated (\"real\") and\nunannotated (“cryptic/novel”) exons, which showed comparable distribution\n(Figure R2.2a).\nRegarding the relationship between DyScore and FAScore, the scatter plot\n(Figure R2.2b) showed a weak correlation coefficient (0.31). Upon close\ninspection, we observed that the high density of events exhibiting a high\nFAScore and low DyScore pattern (as indicated in the red circled area),\nclearly demonstrating that FAScore can capture events with low DyScore that\nmay be missed by filtering based only on DyScore. Furthermore, we plotted the AS\nevents in the positive set, revealing that some known functional AS events had\nlow DyScore values (Figure R2.2c). To further dissect whether the model is\nbiased towards events with high DyScore, we plotted a heat map and found\nthat FAScore and DyScore were not related (Figure R2.2d). These\nobservations further underscore the necessity of using an integrated model for\ncandidate selection, such as FAScore, rather than focus on one or two key\nfeature(s).\nFigure R2.2. Correlation between FAScore and absolute dynamic score across\nalternative splicing events.\na, Density curves show the distribution of FAScore values (ranging from 0.00 to 1.00)\nfor annotated events (yellow) and unannotated events (blue). Both distributions are\nbimodal, with a major peak at low FAScore values (non-functional AS) and a minor\npeak at high FAScore values (functional AS).\nb, Density plot showing the correlation between FAScore and absolute DyScore in\nexon skipping events, with darker regions indicating higher density of events. The red\ncircle highlights a region containing splicing events with high FAScore but low dynamic\nscore. Spearman’s correlation coefficient and corresponding p-value are shown.\nc, Scatter plot showing the correlation between FAScore and absolute DyScore. Each\npoint represents a known functional splicing event. The top five and last five AS events\nare labeled.\nd, Heat map displays normalized feature scores for high DyScore with high FAScore\n(50 were randomly selected) or low FAScore (50 were randomly selected) AS\nevents identified via PU-learning. Each row represents an AS event, and each column\ncorresponds to a predictive feature, with the final FAScore. Values are z-score\nnormalized with blue indicating low values and red indicating high values.\n4. Perhaps more important: the standard in the field of AS to rank candidates\nfor experimental validations is to select evolutionarily conserved events that\nchange their inclusion levels across lineages and/or a time course. Is the\nFAScore approach really much better than this? The authors selected four\nevents with high FASscores and high DyScores (i.e. change in their inclusion\nlevels across lineages and/or the time course). If I understood correctly, the\nmost important contributor to the FAScore is conservation and then lineage\nregulation (which makes it redundant with DyScores). In fact, based on Fig. 2d,\nthe events selected for follow up have a fairly modest FAScore (0.55-0.70) and\nthere are at least dozens of events with much higher FAScores. Why have those\nevents not been prioritized? In other words, the FAScore is perhaps useful, but\nit seems potentially an oversell.\nResponse: As shown above in point 9 in Reviewer 1, we found that the\nconservation and FAScores are weakly correlated (Revised Extended Data Fig.\n3a), FAScores have added value in predicting functional AS events with low\nconservation. Indeed, the functional events with published experimental\nevidence such as Kdm1a-E7, EHMT2-E10 have low sequence conservation\nhowever had high FAScores of 1.00 and 0.98, respectively (Revised Extended\nData Fig. 3b), showcasing the value of FAScore in the selection of AS with low\nsequence conservation, which integrates multiple features, such as\nconservation, dynamic changes and transcript structural composition.\nRegarding the selection of candidate events using FAScore and DyScore, we\nagree with the raised concern and we now updated the selection using species-\nand cell lineage-specific FAScores. Briefly, we selected AS candidates based\non FAScore on a cell lineage in human or mouse (Revised Fig. 3d). For each\nlineage, we showcased one or two events with a high FAScore: erythroid\nlineage, TBC1D23-E15 (0.998) and EPB41L1-E14 (0.995); myeloid lineage,\nKLF6-E3 (0.997) and SSBP3-E6 (0.915); lymphoid lineage, FYN-E7 (0.987),\nwhich was reported to play an important role in T cell signal transduction in\nmice(Kumar, Li et al. 2021). Given their high FAScores in their specific lineage,\nthey were then validated in the experiments.\nRevised Extended Data Fig. 3. Correlation between FAScore and evolutionary\nconservation across alternative splicing events.\na, Density plot showing the correlation between FAScore and evolutionary\nconservationacrossspeciesinESevents,withdarkerregionsindicatinghigherdensity\nof events. The red circle highlights a region containing splicing events with high\nFAScore but low conservation. Spearman correlation coefficient and corresponding p-\nvalue are shown.\nb, Scatter plot showing the correlation between FAScore and evolutionary\nconservation score across species. Each point represents a known functional splicing\nevent.\nRevised Fig. 3d. Rank-ordered distribution of functional AS scores for exon skipping\nevents in five clusters (11, 8, 22, 5, 12). AS events are ranked by descending FAScore\n(0.00 to 1.00). Each point represents an AS event. The number of SE events in each\ncluster is indicated in parentheses. The candidate AS events are highlighted in red.\nThe marginal miniatures in the lower-left corners depict the average FAScore trend of\nevents within the cluster across four hematopoietic contexts: fetal hematopoietic organ\n(FHO), erythroid (Ery), lymphoid (Lym), and myeloid (Mye) lineages. Event sizes (N)\nfor each cluster are annotated below the miniatures.\n5. What gene sets were used as background for the GO enrichment analyses?\nIn the Methods section it is unclear if they have used any background at all. If\nnot, it is possible that the enrichments they obtained were due to lineage-\nspecific gene expression biases, which are known to be a major confounder for\nAS analyses and could explain why they fit the lineages functions well.\nResponse: We appreciate the reviewer’s valuable suggestion. The background\nset for all GO enrichments comprised all genes with GO annotations in the\norganism-specific database (org.*.eg.db). However, we fully acknowledge the\nreviewer's concern that this default background may introduce lineage-specific\nexpression biases, potentially confounding AS-based functional enrichment.\nTo address this, we performed GO enrichment analyses using a new\nbackground gene set consisting of genes that were expressed (TPM >= 0.5) in\nat least one sample in the lineage studied. Lineage-specific functional\nsignatures were still exhibited in the lineage-restricted background (Revised Fig.\n2e). For example, genes harbouring predicted functional AS in erythroid lineage\nconsistently enriched “erythrocyte differentiation”. Genes with predicted\nfunctional AS in enriched in myeloid cell differentiation and neutrophil activation.\nWe revised the ‘Methods’ section (Pg. 29, ln 823-827) to explicitly state:\n“Gene Ontology (GO) enrichment analysis was performed separately for\neach lineage using the R package clusterProfiler95 (v4.2.2). For each lineage,\nthe set of genes harboring functional and uncertain SE events was tested for\nenrichment against a background set comprising all genes that were\nexpressed (TPM >= 0.5) in at least one sample within one lineage to control\nfor cell lineage-specific expression biases”\nRevised Fig. 2e. Bubble plot shows the GO pathways analysis of functional AS-\nharboring genes using cell lineage-specific background sets across FHO development\nand lineage differentiation (including erythroid, lymphoid, and myeloid) in human (Left)\nand mice (Right). The size of each bubble corresponds to the number of genes in the\ninput list that are mapped to a particular GO term. The x-axis represents the\nsignificance of enrichment (p-value < 0.05).\n6. The zebrafish morpholino experiment is much weaker than all the other\nexperiments. In my opinion, it is OK to keep these results, but I would move\nthem after the first results in mouse (i.e. “white blood cell counts remained\nunchanged”) and introduce them in a very specific manner, e.g.: “To assess the\nconservation of this phenotype, we ....”\nResponse: Thanks for the constructive suggestion. We have moved the\nzebrafish results after the results in mouse. And revised the MS as follows (Pg.\n11-12, ln. 318-326).\n“…E15-KO mice exhibited significant reductions in red blood cell (RBC)\nnumbers, hematocrit (HCT) and hemoglobin (HGB) levels compared to the\nWT littermate controls while white blood cell counts remained unchanged\n(Fig. 4e, f). In addition, to assess the conservation of this phenotype, we\ntargeted the long isoform of Tbc1d23 in zebrafish (Tbc1d23-L, including E15)\nby morpholino (MO) treatment. The Tbc1d23-L-MO (E15-MO) treatment\nspecifically reduced the mRNA level of the long isoform while increasing the\nexpression of the short isoform, with no significant effect on total Tbc1d23\nprotein level comparable to the fish treated with the control MO (Ctrl-MO)\n(Fig 4g-i, Extended Data Fig. 7g)…”\n7. In summary, I believe the experimental and mechanistic aspects of this\nmanuscript are much stronger than the FAScore, which is perhaps not needed.\nMy recommendation is to downplay this part, and simply sell it as a way to\nprioritize exons for experimental follow up. As it is now, it has potential flaws\nand its relevance is oversold, which may require substantial work to solve.\nResponse: As mentioned, we updated our model using ‘true’ positive set and\nfocusing on only exon-splicing events and the updated FAScore demonstrated\nimproved precision and accuracy. We believe FAScore is valuable in\nidentification of interesting AS such as TBC1D23-E15 which as no defined\ndomains or motifs and events like Kdm1a-E7, EHMT2-E10 which have low\nsequence conservation.\nReviewer #3:\nHu et al. surveyed the transcriptome of the hematopoietic system in six\nvertebrate species and trained a machine learning model to identify functional\nlineage-specific alternative splicing events. Their approach relies on a Random\nForest model trained on a large dataset of differentially spliced events between\nlineages and leveraging positive unlabeled learning to take advantage of a set\nof 71 validated functional events curated from the literature. The interpretation\nof the model yielded a set of 19 features which were predictive for alternative\nevent functionality. The model was then applied to identify several candidate\nfunctional events in lymphoid, myeloid and erythroid commitment. The authors\ngo on to show that one of their candidates, a vertebrate cassette inclusion event\nin TBC1D23, is important for erythropoiesis using in vivo and ex vivo\nexperiments in the zebrafish and mouse hematopoietic cells. They finally\nprovide mechanistic evidence to demonstrate that the inclusion of the\nalternative TBC1D23 exon 15 in the erythroid lineage decreases binding to its\npartner complex RANBP2/RANGAP1, likely leading the activation of\ntranscription factors involved in erythroid differentiation through lower\nSUMOylation of HDAC1.\nThis study tackles the extent and the role of alternative splicing during\nhematopoiesis. This is an important and poorly studied problem with wide\nranging applications in fundamental and translational immunology. Notably, the\nmanuscript presents three valuable components to the field: 1) a detailed\nsurvey of alternative splicing events across hematopoietic lineages, between\nspecies and between developmental stages; 2) a predictive model for splicing\nevent functionality; 3) the characterization of a novel alternative event in\nTBC1D23 that is key during erythropoiesis.\nThe methodology leading to the developmnt of the FAScore and DyScore is\nreasonable and robust. The methods for the experimental validation of\ncandidate functional events are generally sound, but key experiments and\ndetails lack in some sections, in particular in demonstrating the lineage specific\nof TBC1D23 exon 15 inclusion and differential SUMOylation of HDAC1,\nas described below.\n1. The astronomical number of differential/lineage-specific splicing events\nreported and differentially expressed genes in the TBC1D23-E15KO mouse is\nhard to believe and likely obscures some of the most biologically relevant\nevents. This could be due to lenient filtering criteria and/or specific tools and\nmethods used for transcriptomic analyses. Depending on the version and\nparameters used, SUPPA2 can have a very high false-positive rate. The choice\nof performing scRNA-Seq rather than FACS followed by bulk RNA-Seq to\nidentify DEGs in the TBC1D23-E15KO is puzzling. The authors should at least\ncomment on these points and possibly consider revisiting their filtering criteria.\nResponse: We thank the reviewer for the comments and suggestions. As we\nused 10x Genomics 3’ single-cell RNA-seq followed by NGS, which provides\nonly short reads at the 3’end of the transcripts for differential gene expression\nanalysis; we therefore did not perform differential splicing analysis. We agree\nthat the version and parameters of SUPPA2 can affect the identification of the\nnumber of splicing events, we therefore used SUPPA2 only to identify the type\nof splicing events according to the GTF annotation files (e.g., exon skipping,\nalternative 5’/3’ splicing, etc.). The large numbers of AS events are the overall\ncounts of all exon-skipping exons from multiple lineages and species.\nRegarding the differential gene expression analysis between the WT and\nTBC1D23-E15 KO mouse in the scRNA-seq, the 5,915 unique differentially\nexpressed genes (DEGs) across all 20 cell types were identified using\ncommonly-used single-cell thresholds in Seurat (adj. p<0.05 & |avg_log2FC| ≥\n0.25), with the Pro-E cells showing the highest number of DEGs (3,346, Figure\nR1a), suggesting the KO mouse may have the largest difference in Pro-E cells\nwhen compared to WT mouse.\nTo assess whether the filtering criteria affect our results, we implemented two\nstricter thresholds: |avg_log2FC|≥0.58 (2,492 DEGs, Figure R3.1b) and\n|avg_log2FC|≥1 (802 DEGs, Figure R3.1c). GO enrichment analysis confirmed\nrelatively consistent pathway signatures (Figure R3.1d) across all three\nthresholds (23 pathways, including critical processes like erythrocyte\ndifferentiation, Figure R3.1e). We agree that the |avg_log2FC|≥0.25 cutoff is\nlenient, we therefore used |avg_log2FC|≥1 as the threshold.\nRegarding our choice of scRNA-Seq rather than FACS followed by bulk RNA-\nSeq, the reasons are as follows: 1) we wanted to analyze the cell-types\ninfluenced most by TBC1D23-E15 KO and to compared with FAScore\nprediction, and such information is better provided by scRNA-seq than bulk\nRNA-seq dependent on the sorted cell-types; 2) scRNA-seq enables us to\nperform transcript factor regulatory networks analysis but not bulk RNA-seq\ndata. The analysis provides potential networks of TBC1D23-E15 KO, this may\nprovide insights into the underlying mechanisms TBC1D23-E15 regulation.\nFigure R3.1 Comparison of different threshold values of differentially expressed\ngenes (DEGs) analysis across cell types.\na-c, Bar plots quantifying cell-type-specific DEGs counts under increasing fold-change\nthresholds (adjusted p.value < 0.05): |log FC| ≥ 0.25 (a), |log FC| ≥ 0.58 (b), |log FC|\n2 2 2\n≥ 1.0 (c). Cell types are color-coded, and the number of DEG is labeled on the bar.\nd, Venn diagram illustrating the overlapping of significantly enriched GO pathways of\nDEGs in three thresholds (Biological Processes, adjusted p-value < 0.05).\ne, Dot plot displaying the overlapped pathways in GO-BP enrichment for DEGs\nidentified in (d). The erythroid-related pathway is colored in red.\n2. A lot of the evidence cited to support assumptions is anecdotal. For example,\nthe identification of FYN-E5 and PIK3R1-E1 among high FAScore and DyScore\nevents is used as a demonstration that high FAScore and DyScore can be used\nto find lineage-specific functional events, but how many known lineage-specific\nfunctional events can be found among low FAScore/DyScore events?\nResponse: As mentioned, Reviewer 2 point 1, we agreed that the positive set\nshould only include ‘actual’ functional events with experimental evidence and\ntherefore updated the positive set (Revised Supplementary Table 3). Based\non the updated positive set, we revised a prediction model specifically for exon\nskipping. Moreover, by previously Reviewer 2’s comment (point 4), we revised\nthe selection criteria of candidate AS to rank AS events only based on the\nFAScore in different species without DyScore.\nRegarding the question of how many known functional events have low\nFAScores (i.e., potential false negatives), we included all known functional AS\nevents with published experimental evidence in the training and validation set.\nThe updated model was trained using reliable positive and negative sets\nderived via a PU-learning framework, with 20% of the data reserved for\nvalidation. The trained model achieved a false negative rate of 0.05. As shown\nin Figure R3.2, the true positive set in the validation set have high FAScores.\nAmong these, Mbnl1-E5 predicted to have a score with 0.470 was classified as\nuncertain, although it has been previously reported in the literature to play a\nfunctional role during terminal erythroid differentiation in mice (Cheng, Shi et al.\n2014)\nFigure R3.2 Rank-ordered distribution of FAScore for all AS events in validation set\nfrom PU-learning strategy. AS events are ranked by descending FAScore (0.00 to 1.00).\nThe functional AS events from literature are explicitly labeled.\n3. The experiments in Figure 2E-G must be better described. Currently it is\nimpossible to tell from the main text or figure legend that these experiments use\nmouse progenitor cells that were transduced with exon targeting gRNAs and\nthen sorted by lineage to perform the colony forming assays. Even the methods\ndo not specify that the experiment was performed in mouse cells and what\nmouse strain were used. There should also be a panel showing baseline exon\ninclusion levels by RT-PCR in the control GM, GEMM and BFU-E lineages.\nResponse: Following the reviewer’s suggestion, we have updated the\ncorresponding Results, Figure legends, and ‘Methods’ section.\nFor Results (Pg. 10-11, ln. 270-292):\n“SSBP3-E6 and KLF6-E3 are more frequently included in the myeloid\nlineage (average PSI: 0.27 and 0.92, respectively) compared to the erythroid\n(0.08 and 0.77) and lymphoid (0.10 and 0.75) lineages (Fig. 3e). In contrast,\nTBC1D23-E15 is more likely to be included in the erythroid lineage (average\nPSI: 0.78), while tends to be skipped in the lymphoid (average PSI: 0.32)\nand myeloid (average PSI: 0.39) lineages (Fig. 3e). We first validated that\nall of the candidate exons were all transcribed in the colony-forming unit-\ngranulocyte-macrophage (GM), the colony-forming unit-granulocyte-\nerythrocyte-monocyte-megakaryocyte (GEMM) and the burst-forming unit-\nerythrocyte (BFU-E) colonies derived from colony-forming unit (CFU) assays\n(Extended Data Fig. 6c). To validated that they are functional in lineages\nwhich they have high FAScores, we performed CFU assays on the mouse\nHSPCs in which the candidate exons were deleted using GFP-expressing\nretroviral vector having two sgRNAs targeting the flanking intronic regions of\nthe exons (Fig. 3f). The transduced cells with GFP expression were sorted\nand plated for CFU assays. The CFU assays revealed that the deletion of\nSSBP3-E6 or KLF6-E3 in mouse c-kit⁺ cells resulted in a reduced formation\nof myeloid colonies by 58% and 38% respectively, while the erythroid\ncolonies were of wildtype level (Fig. 3g). Similarly, the CFU assays showed\nthat the deletion of the Tbc1d23-E15 inhibited erythroid colony formation by\n35% without a significant impact on myeloid populations (Fig. 3g). EPB41L1-\nE14 which exhibits a higher inclusion rate in mature myeloid cells (Fig. 3e)\nand the deletion of Epb41l1-E14 markedly enhanced erythropoiesis with an\nincrease of 38% of BFU-E colonies while no significant difference on the\nmyeloid colonies was observed (Fig. 3g). In summary, FAScore accurately\nidentified AS events play important roles in lineage fate commitment.”\nFor Figure legends (Pg. 61-62):\n“Fig.3 | FAscore predicts key AS events that are dynamically regulated\nduring lineage commitment.\ne, Schematic tree view of mean Percent Spliced In (PSI) values for EPB41L1\nexon 14 (EPB41L1-E14), TBC1D23 exon 15 (TBC1D23-E15), KLF6 exon 3\n(KLF6-E3), SSBP3 exon 6 (SSBP3-E6) skipping events derived from RNA-\nseq data of hematopoietic differentiation in human. HSC: Hematopoietic\nstem cell, MPP: Multipotent progenitor, CLP: Common lymphoid progenitor,\nCMP: Common myeloid progenitor, GMP: Granulocyte monocyte progenitor,\nMEP: Megakaryocyte erythrocyte progenitor, EB: Erythroblasts, MK:\nMegakaryocytes, Neu: Neutrophil, Mo: Monocyte, Mp: Macrophage, Mp_act:\nActivated macrophage, Mp_inf: Inflammatory macrophage, B: B cell, CD4:\nCD4+ T cell, CD8: CD8+ T cell, NK: Natural killer cell.\nf. Schematic of CRISPR/Cas9-mediated exon deletion in mouse c-kit+ cells.\nc-kit+ cells isolated from 8-week female Rosa-cas9 knock-in mice were\ntransduced with retroviruses expression two sgRNAs targeting the upstream\nintron (5’ sgRNA) and downstream intron (3’ sgRNA) intronic regions flanking\nthe target exon (indicated by scissors) (Top). PCR was performed using\nexon-KO-F and exon-KO-R primers to confirm exon deletion. Genome PCR\nanalysis of targeted deletion: Epb41l1 exon 4 deletion (Epb41l1ΔE14), Klf6\nexon 3 deletion (Klf6ΔE3), and Ssbp3 exon 6 deletion (Ssbp3ΔE6), and\nTbc1d23 exon 15 deletion (Tbc1d23ΔE15). The red arrows indicate the\ncorrect bands, and all bands were confirmed by Sanger sequencing (Bottom).\ng. Erythroid-myeloid differentiation potential of 5000 control or mutant c-kit+\ncells (Epb41l1ΔE14, Tbc1d23ΔE15, Klf6ΔE3, and Ssbp3ΔE6) from f).\nColony types quantified: BFU-E, colony-forming unit-granulocyte-\nmacrophage (GM), and colony-forming unit-granulocyte-erythrocyte-\nmonocyte-megakaryocyte (GEMM) colonies (n = 3 per group).”\nFor Methods (Pg.33, ln. 964; Pg. 40, ln. 1156-1165; Pg. 34, ln. 976-985):\n“Blood routine, cell isolation and flow cytometry analysis. …HSPCs\nwere further enriched using c-kit/CD117 MicroBeads (Miltenyi Biotec, #130-\n09224)…”\n“Generation of exon-deleted and overexpression mouse HSPCs. The\ntarget exon of Tbc1d23, Ebp41l1, Ssbp3, and Klf6 were deleted by two\nsgRNAs (listed in Supplementary Table 8) were cloned into pMSCV-U6-\nsgRNA-IRES-GFP, which was modified from pMSCV-IRES-GFP (Addgene,\nplasmid #20672). The coding sequence of HDAC1 and HDAC1mut were\ncloned into pMSCV-IRES-GFP for overexpression.\nThe resulting constructs were cotransfected with retroviral packaging\nplasmids (pCL-Eco, pVSV-G) to generate retroviral particles as described\nabove. The viral supernatants were harvested and added into c-kit+ cells of\nCas9+ or WT mice with 8 μg/mL polybrene. At 24 h after transfection, the\ntarget cells were enriched by GFP.”\n“Colony-forming unit assay. Freshly isolated whole BM cells, enriched c-\nkit+ cells or sorted GFP+ cells were plated in methylcellulose culture medium\n(MethoCult™, STEMCELL Technologies) and incubated at 37°C with 5%\nCO . For cells plated on MethoCult™ GF M3334, colonies were counted and\n2\nscored as colony-forming unit erythroid progenitor cells (CFU-E) at 48 h after\nplating. Cells plated on MethoCult™ GF M3436 were evaluated 10–14 days\nafter plating and scored as burst-forming unit erythroid progenitor cells (BFU-\nE). For cells plated on MethoCult™ GF M3434, colonies derived from BFU-\nE and granulocyte-macrophage progenitor cells (CFU-GM, CFU-G and CFU-\nM), multipotent granulocyte, erythroid, macrophage, megakaryocyte\nprogenitor cells (CFU-GEMM) were counted 9–12 days after plating.”\nAdditionally, regarding using RT-PCR to show the baseline exon inclusion levels\nof Tbc1d23, Epb41l1, Klf6, Ssbp3 in GM, GEMM, and BFU-E colonies, we\nwashed down the colonies from plates of MethoCult™ GF M3434 using PBS\nand sorted the cells from these clones by FACS and performed RT-PCR to\ndetermine the inclusion levels of these exons in these cells. In revised\nExtended Data Fig. 6c, we confirmed that the transcription of these AS events\n(with exon inclusion) in the corresponding lineages.\nRevised Extended Data Fig. 6c. Quantification of Percent Spliced In (PSI%) values\nfor Klf6 exon 3, Ssbp3 exon 6, Tbc1d23 exon 15, and Epb41l1 exon 14 was performed\nby RT-PCR analysis in the following cell populations isolated from colony-forming unit\n(CFU) colonies: c-kit+ cells, megakaryocytes (MK), erythroid burst-forming units (BFU-\nE), macrophages, and granulocytes.\n4. The specificity of TBC1D23 exon 15 inclusion is not shown in any of the\nfigures. This is a critical piece of missing data. Evidence should be provided in\nthe form of RT-PCR on several sorted human and/or mouse hematopoietic cell\npopulations.\nResponse: As suggested, we performed RT-PCR using sorted mouse\nhematopoietic cell populations, and the PSIs collated well with the PSIs derived\nfrom bulk RNA-seq analysis of mouse HSPCs (Revised Extended Data Fig.\n7b, c).\nRevised Extended Data Fig. 7b, Box plot of PSI values for Tbc1d23-E15 skipping\nevents derived from RNA-seq data of hematopoietic differentiation in mouse. LT-HSC:\nLong-time hematopoietic stem cell, ST-HSC: Short-time hematopoietic stem cell, MPP:\nMultipotent progenitor, CLP: Common lymphoid progenitor, CMP: Common myeloid\nprogenitor, GMP: Granulocyte monocyte progenitor, MEP: Megakaryocyte erythrocyte\nprogenitor, Ery: Erythroblasts, MK: Megakaryocytes, Mo: Monocyte, Mp: Macrophage,\nB: B cell, CD4: CD4+ T cell, CD8: CD8+ T cell, NK: Natural killer cell.\nRevised Extended Data Fig. 7c, Differential expression of spliced transcripts was\nvalidated in the mouse hematopoietic cells through RT-PCR and corresponding\nagarose gel electrophoresis. Calculate the Percent Splice In (PSI) value of each lane\nwas calculated by dividing the gray value of the longer transcript by the sum of the gray\nvalue of the longer and shorter transcripts.\n5. Does the MO in Fig 3A-B also affect the levels of the TBC1D23 exon skipping\ntranscript? Currently, it is impossible to confirm that the phenotype is due to loss\nof TBC1D23-E15 specifically or loss/decrease of all TCB1D23 expression\naltogether.\nResponse: To address the concern, we measured total TBC1D23 protein\nlevels in MO-treated embryos and found no significant change in total TBC1D23\nprotein expression level (Revised Fig. 4h, i), suggesting E15-MO primarily\naffects the expression of the long isoform but not the exon-skipping short\nisoform.\nRevised Fig. 4h, Immunoblot analysis of total Tbc1d23 protein levels in zebrafish\nembryos following Ctrl-MO and E15-MO treatment.\nRevised Fig. 4i, Images of membranes in (h) were analyzed using ImageJ software to\ncalculate the relative protein level of TBC1D23 (n=3).\n6. Supplementary Fig 5D should include a panel showing TBC1D23 isoform\nratios by RT-PCR using a single primer pair in flanking exons.\nResponse: We appreciate the reviewer’s suggestion. To clarify, the RT-PCR\nanalysis of TBC1D23 isoform ratios (using primers in flanking exons) was\npresented in the original Fig. 3D. However, we acknowledge that this may have\ncaused confusion with the genotyping data. As suggested, we have now\nrestructured the figures as follows: (1) The genotyping PCR results were\nincluded in revised Fig. 4a, (2) the original RT-PCR isoform-ratio data were\nmoved to revised Fig. 4b, and (3) the qPCR results quantifying total Tbc1d23\nmRNA levels were moved to revised Extended Data Fig. 7d.\nRevised Fig. 4a, Schematic representations of the WT and E15-KO alleles, with arrows\nindicating primer binding sites and the expected PCR product sizes labeled (Top) and\ngenotyping of E15-KO mice (Bottom).\nRevised Fig. 4b, RT-PCR validated the specific deletion of the exon 15-containing\nisoform (Tbc1d23-L) in E15- KO mice.\nRevised Fig. 7f. RT-qPCR analysis of Tbc1d23-L (left), and Tbc1d23-S (right) mRNA\nlevels in WT and E15-KO mice.\n7. Nowhere is the difference in HDAC1 SUMOylation between TBC1D23-\nE15KO and WT shown. This is another critical piece of missing data.\nResponse: As suggested, we performed an additional Co-IP analysis of\nSUMOylated HDAC1 using bone marrow cells isolated from both WT and E15-\nKO mice and the results were included in the Revised Fig. 5m-n and MS (Pg.\n14, ln. 401-402):\n“Furthermore, HDAC1 SUMOylation was elevated by 1.5-fold in BMCs from\nE15-KO mice compared to those of WT controls (Fig 5m, n).”\nRevised Fig. 5m, Co-IP analysis of the HDAC1 SUMOylation level in WT and E15-KO\nmice. Bone marrow cells (1×10⁸) were lysed and subjected to IP with anti-HDAC1\nantibody and immunoblot with anti-SUMO1 antibody.\nRevised Fig. 5n, Membrane images in (m) were analyzed using ImageJ software to\nquantify the relative HDAC1 SUMOylation level in E15-KO mice compared to WT mice\n(n=3).\n8. TBC1D23 E15 may alter interaction with partners, but the changes shown in\nco-IP experiments with RANBP2/RANGAP1 are rather small and downstream\nchanges in HDAC1 SUMOylation are not quantified; how can the loss of a\nTBC1D23 isoform be expected to lead to the differential expression of close to\n6,000 genes?\nResponse: As shown above in point 7, we now confirmed the HDAC1\nSUMOylation level was increased by 50% in the bone marrow cells of E15-KO\nmice. To address the concern about interaction the changes of interaction\nbetween TBC1D23 and RANBP2/RANGAP1, we performed an additional Co-\nIP assay using cells expressing only long or short isoform, which were\ngenerated by expressing long or short isoform in the full gene KO cells and the\nexpression level was analyzed by Western blot showing the overall protein\nexpression level was comparable as control cells (Revised Fig. 5c, d).\nQuantitative analysis revealed a 4-fold higher RANBP2 binding by the short\nversus long isoform under these physiological conditions (Revised Fig. 5e, f).\nThis difference was substantially attenuated under overexpression conditions\n(only 1.5-fold difference, p < 0.01), confirming that physiological expression\nlevels are crucial for revealing authentic isoform-specific binding properties.\nRegarding to the transcriptional changes, we repeated the analysis using more\nstringent cut-off (FDR < 0.05, fold-change ≥ 2) and identified 802 differentially\nexpressed genes between E15 KO and WT cells. As HDAC1 is a master\nepigenetic regulator, and it has been reported in the literature that SUMOylation\nof HDAC1 alters its biochemical activity(David, Neptune et al. 2002) and\nenhances its transcriptional repression mechanisms (Yang and Sharrocks\n2004), it is reasonable that the 50% upregulated level of HDAC1 SUMOylation\nleads to scalable changes downstream.\nRevised Fig. 5c, Western blot analysis of TBC1D23 protein expression levels in 293T\ncells under four conditions: sgRNA-control cells, TBC1D23-KO cells reconstituted with\neither empty Flag vector (Flag), Flag-tagged short isoform of TBC1D23 (Flag-\nTBC1D23-S), or Flag-tagged long isoform of TBC1D23 (Flag-TBC1D23-L).\nRepresentative results from three independent experiments are shown.\nRevised Fig. 5d, Statistical analysis of Tbc1d23 expression level. The level was\ndetermined by normalizing the image gray values in (c) (n=3).\nRevised Fig. 5e, Co-immunoprecipitation (co-IP) analysis of the binding capacity of\nTBC1D23-S and TBC1D23-L to the RanBP2 complex. HA-RanGAP1 and EGFP-\nRanBP2-E3 were co-transfected to TBC1D23-KO 293T cells reconstituted with Flag,\nFlag-TBC1D23-S, or Flag-TBC1D23-L. Cell lysates were subjected to\nimmunoprecipitation (IP) with anti-Flag antibody and immunoblot with anti-Flag, anti-\nHA, and anti-EGFP antibodies.\nRevised Fig. 5f, Images of membranes in (e) were analyzed using ImageJ software to\ncalculate the relative intensity of bands for HA-RanGAP1 and EGFP-RanBP2-E3 in IP\nsample (n=3).\n9. Line 209: “classed” should be “classified”.\nResponse: We apologize for the mistake. It was corrected in the revised MS\n“These AS events were then classified as primate-specific (age class 1),\nmammal-specific (class 4) or ancient (present in all vertebrates, class 5,\ndating back 429 million years ago (MYA) (Fig. 3a)\nReferences\nBarbosa-Morais, N. L., M. Irimia, Q. Pan, H. Y. Xiong, S. Gueroussov, L. J. Lee, V.\nSlobodeniuc, C. Kutter, S. Watt, R. Colak, T. Kim, C. M. Misquitta-Ali, M. D. Wilson, P.\nM. Kim, D. T. Odom, B. J. Frey and B. J. Blencowe (2012). \"The evolutionary landscape\nof alternative splicing in vertebrate species.\" Science 338(6114): 1587-1593.\nBekker, J. and J. Davis (2020). \"Learning from positive and unlabeled data: a survey.\"\nMachine Learning 109(4): 719-760.\nCheng, A. W., J. Shi, P. Wong, K. L. Luo, P. Trepman, E. T. Wang, H. Choi, C. B. Burge\nand H. F. Lodish (2014). \"Muscleblind-like 1 (Mbnl1) regulates pre-mRNA alternative\nsplicing during terminal erythropoiesis.\" Blood 124(4): 598-610.\nDavid, G., M. A. Neptune and R. A. DePinho (2002). \"SUMO-1 modification of histone\ndeacetylase 1 (HDAC1) modulates its biological activities.\" J Biol Chem 277(26):\n23658-23663.\nFiszbein, A., L. E. Giono, A. Quaglino, B. G. Berardino, L. Sigaut, C. von Bilderling, I.\nE. Schor, J. H. Enrique Steinberg, M. Rossi, L. I. Pietrasanta, J. J. Caramelo, A.\nSrebrow and A. R. Kornblihtt (2016). \"Alternative Splicing of G9a Regulates Neuronal\nDifferentiation.\" Cell Rep 14(12): 2797-2808.\nHao, Y., R. Colak, J. Teyra, C. Corbi-Verge, A. Ignatchenko, H. Hahne, M. Wilhelm, B.\nKuster, P. Braun, D. Kaida, T. Kislinger and P. M. Kim (2015). \"Semi-supervised\nLearning Predicts Approximately One Third of the Alternative Splicing Isoforms as\nFunctional Proteins.\" Cell Rep 12(2): 183-189.\nHunsicker, M. E., C. V. Kappel, K. A. Selkoe, B. S. Halpern, C. Scarborough, L. Mease\nand A. Amrhein (2016). \"Characterizing driver-response relationships in marine pelagic\necosystems for improved ocean management.\" Ecol Appl 26(3): 651-663.\nKalsotra, A. and T. A. Cooper (2011). \"Functional consequences of developmentally\nregulated alternative splicing.\" Nat Rev Genet 12(10): 715-729.\nKumar, P., L. Li and J. H. Schatz (2021). PI3K/AKT Activation Mediates B-Cell\nTransformation By the\" T\" Splice-Variant of Fyn Kinase, American Society of\nHematology Washington, DC.\nLove, M. I., W. Huber and S. Anders (2014). \"Moderated estimation of fold change and\ndispersion for RNA-seq data with DESeq2.\" Genome Biol 15(12): 550.\nMazin, P. V., P. Khaitovich, M. Cardoso-Moreira and H. Kaessmann (2021). \"Alternative\nsplicing during mammalian organ development.\" Nat Genet 53(6): 925-934.\nMerkin, J., C. Russell, P. Chen and C. B. Burge (2012). \"Evolutionary dynamics of gene\nand isoform regulation in Mammalian tissues.\" Science 338(6114): 1593-1599.\nWang, T., Y. Cao, Q. Zheng, J. Tu, W. Zhou, J. He, J. Zhong, Y. Chen, J. Wang, R. Cai,\nY. Zuo, B. Wei, Q. Fan, J. Yang, Y. Wu, J. Yi, D. Li, M. Liu, C. Wang, A. Zhou, Y. Li, X.\nWu, W. Yang, Y. E. Chin, G. Chen and J. Cheng (2019). \"SENP1-Sirt3 Signaling\nControls Mitochondrial Protein Acetylation and Metabolism.\" Mol Cell 75(4): 823-\n834.e825.\nYang, S. H. and A. D. Sharrocks (2004). \"SUMO promotes HDAC-mediated\ntranscriptional repression.\" Mol Cell 13(4): 611-617.\nZeddies, S., S. B. Jansen, F. di Summa, D. Geerts, J. J. Zwaginga, C. E. van der\nSchoot, M. von Lindern and D. C. Thijssen-Timmer (2014). \"MEIS1 regulates early\nerythroid and megakaryocytic cell fate.\" Haematologica 99(10): 1555-1564.\nZhao, Y., M. C. Li, M. M. Konate, L. Chen, B. Das, C. Karlovich, P. M. Williams, Y. A.\nEvrard, J. H. Doroshow and L. M. McShane (2021). \"TPM, FPKM, or Normalized\nCounts? A Comparative Study of Quantification Measures for the Analysis of RNA-seq\nData from the NCI Patient-Derived Models Repository.\" J Transl Med 19(1): 269.\nZibetti, C., A. Adamo, C. Binda, F. Forneris, E. Toffolo, C. Verpelli, E. Ginelli, A. Mattevi,\nC. Sala and E. Battaglioli (2010). \"Alternative splicing of the histone demethylase\nLSD1/KDM1 contributes to the modulation of neurite morphogenesis in the mammalian\nnervous system.\" J Neurosci 30(7): 2521-2532.\nResponse to Reviewers\nMS Submission ID: NCOMMS-25-31345A\nMS title: The functional landscape of alternative splicing in hematopoietic lineage\ncommitment\nDear reviewers:\nWe thank the reviewers for their valuable comments and suggestions. In response to\nthe comments, we performed additional analysis and revised the manuscript (the\nchanges are shown in red font) and updated the accompanying figures and tables. We\naddressed the concerns point-to-point under each comment. And the changes that we\nintroduced into the revised manuscript (MS) are shown in italics with page (Pg.) and\nline (Ln.) numbers in the revised MS indicated.\nYours sincerely,\nLu Chen, PhD\nReviewer #2\nOverall, the authors have addressed my initial concerns; however, while doing so (and in\nresponse to the other reviewers), I have to say that they have often created more confusion.\nSpecific comments:\n1. Abstract: \"We curated transcriptome data of fetal hematopoietic organ development in\nsix vertebrates and hematopoietic cell differentiation in human and mice, and identified\n84,933 cassette events, originating from 30,255 coding proteins.\" I found this sentence\nconfusing. Checking the main text, it seems it refers to protein-coding genes (as cassette\nexons do not originate from proteins). Thus, I assume this includes all species together (as\nthere are fewer than 30,255 protein-coding genes in human or mouse alone), which is\nmisleading without further context. Also, the exons are likely non-independent (i.e.,\northologous?). This is all too complex for the abstract, in my opinion.\nResponse: We thank the reviewer for pointing this out and great suggestion. We have\nadded an explicit supplementary table (Revised Supplementary Table 2) listing the\nnumber of genes and exons for a specific differentiation lineage and species. The figures\nof 84,933 cassette exons and 30,255 protein-coding genes were obtained by first\nperforming lineage-level deduplication within each species and then summing the counts\nacross different species. We agree that this is too complex, we have revised the sentence\nin the Abstract (Pg. 2, Ln. 35-39) as follows:\n\"We curated transcriptomic data on fetal hematopoietic organ development in six\nvertebrates and hematopoietic cell differentiation in humans and mice. To identify\nfunctional exon-skipping events among thousands of cassette exons in protein-\ncoding genes for a specific differentiation lineage and species, we developed a\nmachine-learning model…\"\nRevised Supplementary Table 2. Number of genes or AS events for a specific\ndifferentiation lineage and species.\n2. In the revised version they claim they focus on 58 experimentally validated cassette\nexons. However, looking at their updated supplementary table (again), there are only 37\nunique exons, as many entries involve the same exon coordinate. And these unique exons\nare actually even fewer, since the human and mouse orthologs are often present in the\ntable. This does not affect the training per se, as this set is only used to define the positive\nand training sets. However, this is a major reporting mistake, to say the least.\nResponse: We have now corrected Supplementary Table 3 by consolidating orthologous\nexons across species and lineages as well as exons with identical genomic coordinates,\nclearly identifying 37 unique cassette exons. Corresponding clarifications have also been\nadded to the main text and Methods section to prevent any potential misunderstanding.\nIn the ‘Results’ section (Pg. 7, Ln. 170-173):\n\"Due to the lack of a verified negative set and the limited number of known\nfunctional AS events with experimental evidence (37 unique exons, with\nredundancies due to cross-species conservation and multi-lineage functionality,\nSupplementary Table 3, See Methods)…\"\nIn the ‘Methods’ section (Pg. 23, Ln. 662-666):\n“A total of 37 uniquely exon-skipping events were experimentally verified as\nfunctional in organ development or hematopoietic differentiation, among 2 AS\nevents observed in two species and 15 observed in multiple lineages or organs\n(Supplementary Table 3), adding up to 58 positive instances. These were used\nto train model totally.”\nRevised Supplementary Table 3. Positive set validated by experiment screened from\nNCBI database.\n3. The use of FAScore throughout the text is confusing and, if I understood it well, a bit\nmisleading. The authors claim that: \"To identify functional exon-skipping events, we\ndeveloped a machine-learning model interrogating 19 features (...), and integrated them\ninto a single prediction score, named FAScore\" (Abstract), and \"Using these features, we\ndeveloped a machine-learning model named FAScore (Functional AS Score) to predict\nfunctional AS events\". If I understood well, there is no single FAScore, but one FAScore for\neach lineage and species, i.e. eight FAScores if I got it correctly. In other words, the text\nshoud be, for instance\" \" To identify functional exon-skipping events *in a given differentiation lineage and species*,\nwe developed a machine-learning model interrogating 19 features (...), and integrated them\nino a single prediction score, named FAScore\" or \"Using these features, we developed\nmachine-learning *models* named FAScores (Functional AS Scores) to predict functional AS\nevents in a given differentiation lineage and species\".\nThe concern is that the term FAScore is used as universal, without this clarification in most\ninstances in the manuscript, which creates confusion. In my opinion, this aspect should be\nexplicitely mentioned and they should refer to a given FAScore as human\nFAScore_Erythroid, FAScore_Lymphoid, etc. And if they pool all FAScores together for all\nlineages (as they often do), this should be made explicit.\nResponse: The reviewer is correct that one FAScore for each lineage and species.\nAccordingly, we have revised the manuscript to eliminate potential confusion throughout\nthe text.\n(Pg. 2, Ln. 37-41)\n“To identify functional exon-skipping events among thousands of cassette exons\nin protein-coding genes for a specific differentiation lineage and species, we\ndeveloped a machine-learning model interrogating 19 features including dynamic\nexpression, protein structure, and evolutionary conservation, and integrated them\ninto a single prediction score, named Functional AS Score (FAScore).”\n(Pg. 6, Ln. 165-170)\n“Using these features, we developed a machine-learning model named FAScore\n(Functional AS Score) to predict functional AS events in a given differentiation\nlineage and species (Fig. 1a). For each AS event, the model assigns a species-\nand lineage-specific FAScore, designated as species_lineage_FAScore (e.g.,\nHuman_erythroid_FAScore, Mouse_erythroid_FAScore).”\n(Pg. 10, Ln. 269-274)\n“Next, we selected candidates with top FAScores from different cell lineages of\nmouse or human for experimental validation (Fig. 3d). For the erythroid lineage,\nTBC1D23-E15 (Human_erythroid_FAScore: 0.998) and EPB41L1-E14\n(Mouse_erythroid_FAScore: 0.995) were chosen. For the myeloid lineage, KLF6-\nE3 (Human_myeloid_FAScore: 0.997) and SSBP3-E6 (Mouse_myeloid_FAScore:\n0.915) were chosen. And for the lymphoid lineage, FYN-E7\n(Mouse_lymphoid_FAScore: 0.987) was selected (Extended Data Fig. 6b).”\nRevised Fig. 3d, Rank-ordered distribution of functional AS scores for ES events in five\nclusters (11, 8, 22, 5). AS events are ranked by descending FAScore (0.00 to 1.00). Each\npoint represents an AS event. The number of ES events in each cluster is indicated in\nparentheses. The candidates AS events are highlighted in red.\nRevised Extended Data Fig. 6b, Rank-ordered distribution of functional AS scores of\nhumans and mice for ES events in cluster 12. AS events are ranked by descending\nFAScore (0.00 to 1.00). Each point represents an AS event. The number of ES events in\neach cluster is indicated in parentheses. The candidate AS event is highlighted in red.\n4. The data points in each boxplot? The FAScores for each age class for each differentation lineage?\nI.e. 4 data points? This would be fine, but it is not clear in the current text and it is not\nexplained in the caption. (Another confusing aspect about this plot: why are there two dots\nin some cases? Also, I would recommend plotting the proportion of funct/non-\nfunct/uncertain for each age category not the other way around.)\nResponse: The data points in each boxplot represent the proportions of AS exons for\ndifferent age classes in each category (FuncAS, NonfuncAS, or Uncertain) in four lineages.\nEach boxplot summarizes four data points, and we have explicitly showed each data point\nusing shape (Figure R1). Furthermore, following the reviewer’s suggestion, we have\nrevised this figure that illustrates the proportions of func/non-func/uncertain AS exons\nwithin each age category (Fig. 3b).\nFigure R1. The percentage of AS exons for different age classes in each category\n(functional, non-functional, or uncertain) in four lineages in human. Each data point\nrepresents the percentage of AS events for an age class within a specific category in one\nlineage.\nRevised Fig. 3b. Proportion of AS exons categorized as functional, non-functional, or\nuncertain across distinct evolutionary age classes of exons in different lineages of human.\n5. Another source of confusion: \"we measured the FAScores for 188,202 annotated and\nunannotated AS events\". Since they replied to my initial comment that they discard\nconstitutive exons, I assume these 188,202 events are for a mix of species or otherwise\nthe numbers do not add up just for humans. Or are they double counting the metrics for\neach FAScore type?\nResponse: As the reviewer correctly understood, 188,202 events indeed represents the\ntotal number of AS events across different species and differentiation lineages used in this\nstudy, We have added a revised Extended Data Fig. 4b and revised supplementary table\n2. (Pg. 7-8, Ln. 199-207):\n“Next, we measured the FAScores for annotated and unannotated AS events\nfrom fetal hematopoietic organ development and hematopoietic cell differentiation\nper species (Supplementary Table 2), that were not used in the model training\nderiving a continuous score for AS, ranging from 0 to 1, with 1 as most likely to\nbe functional. By fitting a three-component GMM, we defined the following classes\nby probability density intersection method: ‘functional AS’ (FuncAS), ‘non-\nfunctional AS’ (NonfuncAS), and ‘uncertain AS’ (Uncertain) (Extended Data Fig.\n4a). In total, we predicted 30,615 functional AS (16.3%), 53,454 non-functional\n(55.3%), and 104,133 uncertain (28.4%) (Extended Data Fig. 4b,\nSupplementary Table 2). ”\nRevised Extended Data Fig. 4b. Stacked bar plot showing the functional classification of\nAS events across all datasets. The y-axis indicates the proportion of events. Colors\ncorrespond to categories as functional, non-functional, or uncertain. The number in plot\nindicates the number of AS events per category.\nRevised Supplementary Table 2. Number of genes or AS events for a specific\ndifferentiation lineage and species.\n6. DyScore: this is now more intuitively explained in the first section of Results. However,\nas above, since the DyScore is specific to a given lineage differentiation (Erythroid,\nLymphoid or Myeloid), it should be clear which lineage they are referring to. For example,\nI assume the cases in Fig. 1 relate only to Erythroid differentiation. If so, this should be\nmore clearly state, and call it e.g. DyScore_Erythroid to avoid confusion. If my\nunderstanding it is not correct, then the authors need to clarify this well to avoid\nmisunderstandings.\nResponse: We have added relevant lineage information to the description of the dynamic\nfeatures, particularly for the cases presented in Fig. 1. The specific modifications are as\nfollows:\n(Pg. 5, Ln.138-139)\n“We quantified six features of the dynamic gene expression of MEIS1 in human\nerythroid lineage.”\nRevised Fig1b-e.\nb, Schematic representation of dynamic profiling for gene expression and alternative\nsplicing (AS). A linear fit (blue dashed line) and a Loess curve (green solid line) are shown,\nwith the shaded area representing the 95% confidence interval. Derived features include:\nrange of genes (Gene.range) or AS (AS.range) to assess variation breadth; correlation\ncoefficient (Gene.cor, AS.cor) and its p-value (Gene.cor.p, AS.cor.p) to quantify association\nwith temporal or differentiation trajectories; linear regression slope (Gene.slope, AS.slope)\nand its p-value (Gene.slope.p, AS.slope.p) to estimate rate and significance of change;\nand tau (Gene.tau, AS.tau) to evaluate specificity.\nc-e, Expression and splicing patterns of MEIS1 and its two AS events during human\nerythroid differentiation. (c) Expression values (transformed TPM) of MEIS1 (Left) and the\nvalues of dynamic features (Right). Percent spliced in (PSI) values (Left) and dynamic\nfeature scores (Right) for the exon 12 skipping event in transcript MEIS1-201 (d) and the\nexon 1 skipping event in transcript MEIS1-202 (e). A linear fit (blue dashed line) and a\nLoess curve (green solid line) are shown, with the shaded area representing the 95%\nconfidence interval.\n7. Perhaps more relevant, since the DyScore is integrated into the FAScore and not used\nanymore beyond Fig. 1, why is it even needed?\nResponse: We reckon DyScore is a easier for readers to understand, providing a single\nvalue to measure the level of dynamic change. These features of DyScore is indeed\nintegrated into the overall FAScore, we calculated the DyScore primarily to quantify the\ndynamicity of genes/AS events. This quantification facilitates the dynamic categorization\nand comparison of genes/AS events, thereby helping to elucidate the importance of\ndynamic changes for functional prediction. Besides, DyScore can characterize the dynamic\ndegree at both the gene and AS levels, providing a useful reference even for those focusing\nsolely on gene-level functionality.\nMinor:\n8. Fig 2c: the FAScore_Lymphoid seems to do quite poorly. Is it the case?\nResponse: We thank the reviewer for this thoughtful observation. Our previous figure\nshowed the difference in splicing site conservation between functional and non-functional\nAS events across various lineages. We update the figure by showing the original\nconservation score from functional and non-functional AS events. As shown in the figure\nbelow, consistent with all other lineages studied, functional AS events in the lymphoid\nlineage also exhibit significantly higher conservation around their splicing sites compared\nto non-functional AS events (Revised Fig. 2c). The lesser difference in the lymphoid\nlineage is due to the conservation level for non-functional AS events in this lineage is\nnotably higher. We have revised the Fig. 2c to improve the clarity.\nRevised Fig. 2c. Comparison of conservation scores (PhastCons) surrounding splice sites\nbetween functional and non-functional AS of annotated type in different lineages. The 3’SS\nand 5’SS refer to the 3’ and 5’ splice sites, respectively. The x-axis indicates the position\n(bp) relative to the splice site, with “0” representing the splice site, while the regions to the\nleft and right are intronic (100 bp) and exonic (50 bp). The y-axis shows the average\nPhastCons scores for each site of functional and non-functional AS.\n9. Extended data Fig. 6: as far as I understood, Mfuzz was designed to cluster continuous\n(temporal) data series, not unrelated classes as those here. If so, this implies that the order\nin which the independent classes are provided could affect the results, which would not\nmake sense.\nResponse: The Mfuzz package, which based on fuzzy c-means clustering, was designed\nfor time-series or continuous gradient data, and its algorithm considers the order of data\npoints. In our study, we used the core fuzzy c-means clustering algorithm (Hathaway and\nBezdek 1986). Notably, fuzzy c-means algorithm is fundamentally a general-purpose\nmachine learning method that does not mandatory require time- or order-based continuity\nbetween data points, and can be used for category data. The essence of this algorithm is\nto cluster data points based on their \"distance\" in a multidimensional space, which is\nindependent of whether the dimensions represent time points. In our analysis, we treated\ndifferent independent classes (FHO, erythroid, lymphoid and myeloid) as distinct\ndimensions in the multidimensional space, and the clustering aimed to identify groups of\nAS events with similar functional patterns across these dimensions.\nTo ensure the robustness and objectivity of our results, we altered the input order of the\ncategories and re-ran the clustering. We found that the AS composition of each cluster and\ntheir function trend patterns remained highly consistent (Figure R2). This indicates that\nour clustering results are robust and their biological interpretation is not dependent on the\ninput order of the categories.\nTo better match the intent of our analysis, we have replaced the original continuous line\nchart with a categorical box plot adding a red line which connects the median values of\neach category to clearly outline the overall trend across the different groups (Revised\nExtended Data Fig. 6a). This change was made to more accurately represent the\ncategorical nature of our data and to avoid confusion of continuity of different lineage.\nFigure R2. Comparison of fuzzy c-means clustering profiles under different category input\norders. Shown are the patterns of four representative clusters (5, 8, 11, 22). Data point\nlabels correspond to the number of AS events per cluster.\nRevised Extended Data Fig. 6a. The distribution of the normalized average FAScore of\nancient AS events (age class 5) for four hematopoietic lineages in 24 clusters from fuzzy\nc-means clustering. Plot numbers indicate the number of AS events assigned to each of\nthese clusters. A red line connects the median values of the lineages, delineating the\noverall function trend across the four lineages for that cluster.\nReviewer #3\nHu et al. provided a substantially revised version of their manuscript. Although none of the\ninitial findings or claims had to be revisited, there were important pieces of evidence and\ntechnical details missing from the original submission. The authors appropriately clarified\nor modified the methodology the address the reviewers' comments. They also provided\nkey results to support their initial claims. Most notably, they explained better how they came\nto such a large amount of differentially expressed/spliced genes and re-ajusted their\nanalysis parameters to highlight the most striking changes, they revisited their prediction\nmodel for exon skipping, they provided baseline quantifications of TBC1D23 isoforms\nacross all relevant hematopoietic cell types, they added missing controls in many\nexperiments, and they substantiated their claims linking the change in TBC1D23 isoform\nratios to HDAC1 SUMOylation through the RanGAP1-RanBP2 complex with standard\nbiochemical and biophysical assays.\n1.There are the multiple typos and incorrect statements (e.g. a mention of \"30,255 coding\nproteins\" in the abstract - but what are \"coding proteins\"?) that populate the text and that\ncan easily be fixed.\nWe believe that the revised manuscript delivers the main message of the authors - that\nthey were able to draft a summary of functional alternative splicing events during\nhematopoiesis - much more successfully than the original manuscript.(Remarks on code\navailability)\nResponse: We appreciate the reviewer’s suggestion and apologize for the mistake. We\nhave now thoroughly proofread the entire manuscript to correct all typos and inaccurate\nstatements. Specifically, we have revised the sentence in the Abstract (Pg. 2, Ln. 35-39)\nas follows:\n\" We curated transcriptomic data on fetal hematopoietic organ development in six\nvertebrates and hematopoietic cell differentiation in humans and mice. To identify\nfunctional exon-skipping events among thousands of cassette exons in protein-\ncoding genes for a specific differentiation lineage and species, we developed a\nmachine-learning model…\"\nReferences\nHathaway, R. J. and J. C. Bezdek (1986). \"Local convergence of the fuzzy c-Means algorithms.\" Pattern\nRecognition 19(6): 477-480."}]

Decision letter