Splice isoform-perturbation coupled to single cell transcriptome profiling reveals functions of microexons in neurogenesis and autism-linked pathways

Domain: neuroscience
Published: 2026-02-03
Editorial outcome: major_revision

Abstract
A major goal of biomedical research is to assign functions to the myriad alternative RNA and protein isoforms. This challenge is particularly relevant to the mammalian nervous system, which produces complex repertoires of alternative splicing events. Here, we describe CHyMErA-seq, a platform that couples systematic deletion of exons to a single cell transcriptomics read-out, and apply this method to investigate a critical program of brain-specific microexons. Perturbation of microexons during neurogenesis reveals convergent roles in the temporal regulation of gene expression programs that direct signaling pathways and morphogenesis. We further observe microexons, including those in the Bin1, Clasp1, Gfra1, Med23, Ptprf and Ralgapb genes, that are required for the correct timing of autism-linked gene expression. Collectively, we describe a flexible system for isoform-resolution perturbation at a single cell level, together with insights into the roles of microexons in the developmental timing of neurogenesis transcriptomic signatures linked to brain disorders.

Introduction
Alternative splicing represents a critical layer of gene regulation that greatly diversifies the expression of transcript and protein isoforms from a limited repertoire of genes. Transcripts from over 95% of human multi-exonic genes undergo alternative splicing1,2, and more than half of these genes encode alternative splicing events that are subject to pronounced condition (e.g., cell/tissue type) dependent regulation3. These events are typically controlled by combinations of RNA binding proteins that interact with cis-elements in precursor (pre-) RNA to promote or repress the proximal assembly of small nuclear ribonucleoprotein particles and additional factors that form active splicing complexes (spliceosomes)4–6. While significant progress has been made towards understanding the regulatory networks that control alternative splicing, a major challenge confronting the post-genomic era is to determine the functions of the large repertoires of uncharacterized isoforms generated by alternative splicing and other post-transcriptional processes that are associated with normal and disease physiology.

Functional characterization of an individual alternative splicing event can require years of effort and typically involves tailored assays that lack scalability. To begin to address this challenge, recent studies have described the development and application of CRISPR-Cas editing strategies in screen formats to assess the impact of exon deletion on cell fitness7–9. For example, ‘paired guide RNAs for alternative exon removal’ (pgFARM) has been used to assess the function of premature termination codon (PTC)-exons that elicit nonsense mediated mRNA decay on cell growth8. We previously described ‘CasHYbrid for Multiplexed Editing and scReening Applications’ (CHyMErA) and applied this system to assay the impact of deleting thousands of alternative exons in two different cell lines on cell fitness7. CHyMErA employs simultaneous expression of Cas9 and Cas12a nucleases, together with hybrid guide (hg) RNAs, which are fusions of Cas9 and Cas12a guides expressed from a single U6 promoter and processed into pairs of individual Cas9 and Cas12a gRNAs that direct genomic fragment deletion. CHyMErA has recently been expanded to a larger screen revealing that ‘fitness exons’ are concentrated in highly expressed essential genes, including important examples such as an alternative exon in the TAF5 subunit that regulates assembly of the TFIID transcription initiation complex and global gene expression9. Also recently, gRNA-directed RNA targeting of deactivated (d)Cas13d (dCasRx) fused to splicing effectors has been used to direct the specific and efficient activation or repression of endogenous exons10,11, although this strategy has not yet been employed on a large scale.

To date, splice isoform-resolution functional genomic assays have not been applied in the context of developmental biology, nor on a single cell level, yet these applications are critical for determining the functions of alternative splicing events linked to cell fate decisions, differentiation programs, as well as cell and tissue maintenance. They are also needed to understand the roles of splice variants and other isoforms in pathological contexts, for example, the numerous forms of cancers and brain disorders in which alternative splicing is dysregulated and in which individual alternative exons are known to contribute to pathophysiology5,6,12. These challenges are particularly relevant to a program of evolutionarily conserved, 3-27 nucleotide-long ‘microexons’, which are predominantly included in the nervous system13,14. Neuronal microexons are enriched in genes that have important roles in nervous system biology and disorders, and they are often misregulated in the brains of individuals with idiopathic autism spectrum disorder (ASD)13. Mice haploinsufficient for the ‘master’ regulator of brain-specific microexons, the Serine/Arginine-repetitive matrix-4 (Srrm4) protein (previously known as nSR10015), which promotes the splicing of the majority of alternative microexons during neurogenesis, exhibit multiple ASD-like phenotypes, including altered social behaviors, hypersensitivity to environmental stimuli, and altered synaptic transmission16,17. Moreover, focused characterization of individual brain-specific microexons has revealed important roles in the regulation of neuronal-specific gene expression18,19, translation20, mRNA deadenylation21, neuritogenesis16,22,23, and synaptic biology24,25. Mice deficient of some of these microexons exhibit neurodevelopmental and behavioral deficits, including ASD-like phenotypes20,21. However, despite the importance of neuronal microexons revealed by these studies, hundreds of these alternative splicing events await characterization.

In this study, we have developed CHyMErA-seq, an enhanced CHyMErA system that enables scalable isoform-resolution perturbation screens coupled to transcriptomic read-outs. We apply CHyMErA-seq to determine roles of neuronal microexons in genes with critical functions in neurodevelopmental processes, and that are implicated in brain disorders. Through applying this system in the context of neurogenesis, we have identified a subset of microexons that convergently control signaling pathways and developmental timing of transcriptomic signatures responsible for critical stages of neuronal maturation and function. Overall, our results introduce CHyMErA-seq as versatile platform for isoform-resolution functional genomics in the contexts of developmental and single cell biology.

Results
A platform for exon-resolution functional genomics at a single cell level
To establish a system enabling the profiling of single cell transcriptomic phenotypes as a consequence of exon perturbation, we coupled CHyMErA-directed deletion of microexons to split-pool single nuclear (sn)RNA-sequencing using the sci-RNA-seq3 protocol26. Specifically, we adapted the ‘CRISPR droplet sequencing’ (CROP)-seq lentiviral vector27 by integrating hgRNA (refer to Introduction) expression cassettes within the vector’s 3ʹ Long Terminal Repeat (LTR). Following genomic integration and LTR duplication, the modified CROP-seq vector expresses hgRNAs from an upstream RNA pol III-transcribed cassette that directs genome editing. The second (duplicated) hgRNA cassette is transcribed by RNA pol II within the 3´ untranslated region (UTR) to generate polyadenylated (polyA) transcripts. When polyA+ scRNA-seq is performed, both cellular and hgRNA transcripts are thus simultaneously profiled, allowing detected transcriptomic changes to be linked to specific microexons targeted for deletion (Fig. 1a).

To initially test and optimize the efficiency of this system, CROP-seq vectors expressing six different hgRNAs targeting three microexons were introduced into a set of CGR8 murine embryonic stem cell (mESC) lines engineered to stably express Cas9 in addition to a series of Lachnospiraceae bacterium (Lb)Cas12a or Acidaminococcus sp. (As)Cas12a constructs. Each set of LbCas12a and AsCas12a constructs incorporated different configurations of nuclear localization signals (NLSs), which were evaluated for nuclear localization efficiency (Yuxi Xiao and J. M., unpublished data), and mutations shown previously to increase the activity of Cas12a28 (Figs. 1b, c; and Supplementary Figs. 1a and 1b). We observe that mutationally-enhanced (en)AsCas12a29 fused to two tandem, carboxy-terminal c-Myc NLSs, expressed together with Cas9, affords the highest overall efficiency of editing compared to the original CHyMErA system, achieving on average a 16% higher editing efficiency (Fig. 1b). Moreover, proportions of edited cells remain constant during neuronal differentiation of the mESC lines, indicating that the microexon deletions, and more generally the introduction of genomic double-stranded breaks, did not significantly affect cell viability (Fig. 1d). The resulting enhanced CHyMErA system was therefore used in downstream experiments to target a larger set of neuronal microexons and analyze the effects of their deletion on differentiating neural cell transcriptomes using snRNA-seq.

A microexon-targeting hybrid guide RNA library
Microexons targeted for analysis by CHyMErA-seq were selected on the basis of their: (1) sequence conservation between human and mouse, (2) exhibiting consistent alternative splicing patterns during in vitro and in vivo neurogenesis30,31, (3) location in genes that have important functions in nervous system biology (Supplementary Fig. 1c), and (4) the availability of multiple gRNA targeting sites in the flanking introns (Fig. 1e, f). All of the targeted microexons are frame-preserving and expected to insert 3-9 amino acids into the coding sequence. Cas9 and Cas12a gRNAs were designed such that at least 100 base pairs (bp) of upstream and downstream sequence flanking the microexon were targeted for deletion, to ensure removal of both the coding sequence and critical splicing signals, while maintaining a distance of at least 150 bp (median distance = 2023 bp) from the upstream and downstream neighboring exons to avoid disruption of their splicing (Fig. 1g). We scored potential Cas9 gRNAs using CRISPOR32 and Cas12a guides using CHyMErA-NET7 (Supplementary Fig. 1d). To ensure robust phenotype detection within a single experiment, for each microexon of interest we selected nine distinct pairs of gRNAs with the highest overall combined scores for targeting. For comparison and as positive controls for possible phenotypic change, we included Cas9 and Cas12a gRNAs that target constitutive exons in the same set of genes. As additional positive controls, we included gRNAs targeting constitutive exons in Srrm4 and its paralog, Srrm3. As negative controls, we included two gRNA pairs designed to cut within intronic sites that are at least 300 bp from the target exon 3’ and 5’ splice sites, and which are therefore not expected to disrupt microexon splicing. Additional negative controls included 13 non-targeting gRNAs, and 13 gRNAs designed to cut within intergenic regions. Collectively, the resulting library consisted of 500 gRNAs targeting 37 microexons in 32 genes, with four genes (Clasp1, Ptk2, Ptprf, and Vav2) containing more than one targeted microexon. We transduced CGR8 mESCs expressing the enhanced CHyMErA system with the lentiviral hgRNA library at an approximate multiplicity of infection of 0.2, such that each cell expressed on average a single hgRNA. To confirm library representation and uniformity, we sequenced the integrated hgRNA cassette from transduced cells after selection. This revealed 100% library representation, with 98.8% of guide sequences having read counts falling within 10-fold of the mean, indicating high uniformity of library coverage (Supplementary Fig. 1e).

Detection of microexon deletion-associated neuronal transcriptomic signatures
To investigate the functions of the targeted microexons during neurogenesis, mESCs expressing the CHyMErA-seq system were differentiated into glutamatergic neurons, and at day in vitro 7 (DIV7) scRNA-seq libraries were generated (Fig. 1h). ~21,000 quality-filtered nuclei expressed targeting or control gRNAs, representing 36 of the 37 microexons targeted by the hgRNA library (Methods). Each pair of microexon-targeting gRNAs was represented by a median of 44 cells, and each microexon target was represented by a median of 492 cells (Supplementary Figs. 2a, and 2b). We detected diverging differentiation trajectories in the transcriptomic profiles and categorized cells into distinct lineages using CellRank233, an analytical framework that computes cell-fate transition probabilities and terminal states from gene expression profiles in scRNA-seq data. Mapping the single-nuclei transcriptomes from an in vivo snRNA-seq dataset from developing mouse embryos26 onto our screen revealed that most cells represent the neuronal lineage, while smaller subsets of cell populations represent early mesodermal and endodermal lineages (Fig. 2a). We clustered cells and analysed marker gene expression within these trajectories to reveal the presence of expected neuroepithelial, radial glial, and neuronal clusters within the neuronal lineage (Fig. 2b, c Supplementary Figs. 3a and 3b), with cells expressing microexon targeting guides for 24/36 detected microexon targets represented by at least 25 cells in each cluster within the neuronal lineage (Supplementary Fig. 3c).

To identify transcriptomic signatures associated with microexon deletion, we analyzed differential gene expression (DGE) profiles of cells representing subclusters within the neuronal lineage (Supplementary Fig. 4a). Specifically, we used a pseudobulk approach34 to calculate pairwise DGE values between subclusters of cells expressing gRNAs directing microexon deletion and non-targeting control gRNAs. Since we expected to only observe transcriptomic signatures in sub-populations of cells with effective gRNA editing, we developed an analytical pipeline to discriminate profiles that likely represent edited (‘knockout’) from unedited (‘non-perturbed’) cells based on these signatures (Methods) (Supplementary Fig. 4b). From this analysis, ten of the analyzed microexons are associated with transcriptomic phenotypes, represented by at least three different pairs of gRNAs targeting the same microexon (Methods) (Supplementary Fig. 4c). Supporting the reliability of these results, benchmarking of our analysis pipeline against a published single-cell CRISPR screen dataset35 reveal that 12/15 (80%) of perturbations called by our method have significant levels of reduction in expression of the target gene (P<0.05, Bonferroni-corrected Wilcoxon rank-sum test).

To further assess specificity of the effects of deleting targeted microexons, we quantified the degree of correlation between the transcriptomic signatures detected in cells expressing distinct gRNAs targeting the same microexon and compared this to correlations between these gRNAs and either intergenic negative controls, or non-targeting controls. Importantly, cells with guides targeting the same gene exhibited significantly higher gene expression correlation compared to both non-targeting (P=7.50x10-12, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) and intergenic-targeting controls (P=5.85x10-9, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) (Supplementary Fig. 4d). This degree of correlation is to a large extent independent of the number of cells recovered for each guide (Supplementary Fig. 4e). We also observe that the expression levels of microexon-harboring genes are not significantly altered upon microexon deletion, indicating that the observed phenotypes are due to microexon deletion rather than loss of gene expression (Supplementary Fig. 4f). To assess the specificity of microexon-deletion phenotypes, we performed pairwise comparisons of microexon deletion-associated transcriptomic signatures between cells expressing microexon-targeting hgRNAs or cells expressing guides targeting constitutive exons in the same gene, with cells expressing non-targeting guides (Methods) (Supplementary Fig. 4g). This analysis reveals that transcriptomic signatures derived from microexon-deletions are distinct from those detected in response to guide-directed perturbation of the corresponding genes. However, we do detect a similar rate of transcriptomic perturbations for gene (i.e., constitutive exon) knockout gRNAs, although these changes involve a different subset of genes from those showing microexon deletion-dependent perturbations (Supplementary Fig. 4h). Moreover, gRNAs targeting Srrm3, and simultaneously Srrm3-Srrm4, also result in a transcriptomic phenotype, consistent with the role of these factors in the widespread regulation of the neuronal microexon program. Intergenic control guides did not generate a phenotype, and negative control intronic deletion gRNAs yielded only a single significant transcriptomic profile change, which may relate to the partial retention of the target intron (in the Fbxo25 gene) during neurogenesis. Overall, these results support the efficacy of CHyMErA-seq for the discovery of specific transcriptomic phenotypes as a consequence of deletion of neuronal microexons during neurogenesis.

Microexon deletions impact transcriptomic signatures associated with neurogenesis
To further investigate the functional consequences of microexon disruption, we performed higher-resolution cell clustering and observed enrichments of perturbed cells at different stages of neurogenesis (Figs. 2c and 2d). This revealed that deletion of a subset of microexons results in altered transcriptomic profiles enriched in subclusters of cells corresponding to late neuroepithelial stages (Bin1, Clasp1(A), Gfra1, Med23, Ptprf(A), Ralgapb), radial glia (Ptprf), and subsequent stages of neuronal differentiation (Bin1) (Fig. 2d). We observe that differentially expressed genes associated with these microexon deletions are enriched in gene sets involved in nervous system development and neuron projection guidance (Fig. 2e). Notably, several of the microexons with transcriptomic perturbations in this study (Bin1, Clasp1, Gfra1, Ptk2 and Ptprf) are located within a network of microexon-containing genes enriched in functions relating to neuron projection development (Fig. 2f). These microexon-containing genes include Ptk2, for which we observe an overall transcriptomic signature associated with deletion of its microexons. Ptk2 encodes focal adhesion kinase (FAK), a non-receptor protein tyrosine kinase that has been shown to mediate focal adhesion dynamics that regulate cell migration and survival in neurons and other cell types36–38. The Ptk2 microexons (denoted as ‘A’, ‘B’ and ‘C’), which display coordinated splicing patterns in vivo, are situated proximal to phosphorylated regions of FAK and have been shown previously to promote FAK phosphorylation39–41. These results imply that microexons have important roles in neuronal differentiation through modulating key proteins involved in morphological changes associated with neurogenesis.

To enable the evaluation of transcriptomic signatures in the context of neurogenesis, we utilized a CellRank2-derived score for each gene according to the correlation of its expression with progression along the neuronal lineage33,42 (Methods). Gene Set Enrichment Analysis (GSEA) reveals that high-scoring genes are enriched in GO terms associated with synapse maturation and related neuronal processes, whereas low scoring genes are enriched in proliferation-related terms, such as chromosome segregation and catabolism (Supplementary Fig. 4i). To assess how microexon deletion effects neurogenesis, we compared the distributions of scores for up- and down-regulated genes for each perturbation, or in the case of Ptk2, pooled scores for genes affected across the three microexon perturbations, due to the relatively low cell numbers and sparse data recovered for each individual Ptk2 perturbation. Notably, compared to background genes (i.e., non-differentially expressed, FDR≥0.05), we observe that upregulated genes associated with altered transcriptomic signatures are significantly enriched for higher CellRank2 scores for microexon deletions in the Bin1, Clasp1, Gfra1, Med23, Ptprf and Ralgapb genes (P<0.05, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) (Fig. 2g). This observation suggests that deletion of these microexons results in premature activation of gene expression programs associated with neurogenesis.

The Gfra1 microexon regulates neuron projection development
We next further explored the functional consequences of deletion of specific microexons that result in the upregulation of neurogenesis gene expression patterns. We initially focused on a microexon in the gene encoding the Glial derived neurotrophic family receptor alpha-1 (GFRA1), which has a prominent role in controlling neuronal differentiation, maturation and signaling. GFRA1 is a glycosylphosphatidylinositol (GPI)-linked surface co-receptor for glial derived neurotrophic factor (GDNF) and known to mediate various signaling pathways through interactions with different surface-exposed co-receptors, including the ‘Rearranged during transfection’ (RET) receptor tyrosine kinase and Neural cell adhesion molecule (NCAM)43,44. GFRA1 contains a single microexon located within a loop region between its first and second GDNF-binding domains45. This region has been associated with reducing GDNF-mediated signaling46, an event that is critical for neurite formation43,44, thus highlighting the importance of characterizing the function of the Gfra1 microexon during neurogenesis.

Analysis of the Gfra1 microexon deletion-associated transcriptomic signature reveals a premature upregulation of programs of genes involved in neurogenesis, including those associated with neuronal projection development (Fig. 3a and Supplementary Fig. 5a). To validate and extend these findings, we generated independent mESC lines harboring homozygous deletions of the microexon (Supplementary Fig. 5b) and differentiated these cells into neurons. From an analysis of bulk RNA-seq from DIV7 neurons generated from these cell lines, we observe a significant positive correlation of log-fold changes in expression when comparing measurements derived from bulk RNA-seq and those derived from the CHyMErA-seq data (R=0.14, p=0.0057, Spearman correlation) (Supplementary Fig. 5c). Furthermore, we observe a significant overlap between the sets of Gene Ontology (GO) terms enriched for upregulated genes in response to Gfra1 microexon deletion in both bulk RNA-seq and CHyMErA-seq data (odds ratio = 20.4, P<2.2x10-16, Fisher’s exact test) (Fig. 3b). We observe that most of the overlapping GO terms are related to processes involved in cellular morphogenesis and neuron differentiation. Additionally, when ranking differentially expressed genes in the bulk RNA-seq with the same CellRank2-derived neuron differentiation score as was used to analyse the CHyMErA-seq signatures, we observe a significant increase in score ranking for upregulated genes (P=2.9x10-4, two-sided Wilcoxon rank sum test) (Fig. 3c). Therefore, to investigate whether the Gfra1 microexon deletion affects neuritogenesis, we visualized and measured neurites by immunofluorescence microscopy using anti-MAP2 antibody to stain neuronal microtubules at DIV3 and DIV7, which correspond to time points representing immature and maturing neurons, respectively (Fig. 3c, and Supplementary Fig. 5d). Notably, we observe formation of longer neurites in DIV7 Gfra1ΔMIC/ΔMIC neurons, compared to wildtype neurons (Fig. 3d; P=0.0056, two-tailed unpaired t-test). Collectively, these results support a role for the Gfra1 microexon in suppressing a transcriptomic program to prevent premature activation of neuritogenesis.

Functions of a neuronal microexon in neurogenesis via modulation of FAK signaling
We next performed a focused follow-up investigation of a previously uncharacterized microexon in the Cytoplasmic linker associated protein 1 (Clasp1) gene, the deletion of which resulted in a pronounced transcriptomic signature in the CHyMErA-seq screen (Fig. 2e). CLAPS1 is a microtubule-associated protein that stabilizes the plus ends of growing microtubules to facilitate directional migration47. In humans, CLASP1 has been genetically linked to ASD through the presence of missense mutations along its coding sequence48. To confirm and extend the results from the CHyMErA-seq screen, we compared altered transcriptomic profiles detected following deletion of the Clasp1 microexon with those detected following direct manipulation of the exon in transcripts using a recently described11 dCasRx-artificial splicing factor fusion protein, dCasRx-RBM25. We have demonstrated that direct gRNA-targeting of dCasRx-RBM25 to downstream intronic sequences proximal to alternative exons lead to their efficient inclusion, whereas targeting dCasRx-RBM25 to exons inhibit their splicing. mESCs stably expressing dCasRx-RBM25 were transduced with gRNAs directing activation or repression of the Clasp1 microexon in mESCs (DIV-8), and during differentiation (Fig. 4a). Bulk polyA+ RNA-seq data was generated and analyzed for the mESC and DIV7 samples (Supplementary Fig. 6a).

Interestingly, we observe premature upregulation of neuronal genes upon forced skipping of the microexon in mESCs compared to when the microexon is subject to forced inclusion (P=6.7×10-6, two-sided Wilcoxon rank sum test) (Fig. 4b). This DGE signature becomes more pronounced at DIV7 (P=9.3×10-27, two-sided Wilcoxon rank sum test) (Fig. 4c), and significantly correlates with the DGE signature detected following deletion of the Clasp1 microexon in the CHyMErA-seq screen (Fig. 4d; R =0.55, P=1.3×10-9, Fig. 4e; R =0.27, P=0.0049. Pearson correlation). We further observe a significant overlap of GO terms enriched among genes that are upregulated in the CHyMErA-seq gene signature and in mESCs (odds ratio = 29.3, P<2.2×10−16, Fisher’s exact test), and DIV7 neurons (odds ratio = 52.1, P<2.2×10−16, Fisher’s exact test) (Fig. 4f). These GO terms are associated with functions related to the regulation of nervous system development, morphogenesis and cell-signaling.

These findings led us to hypothesize that the Clasp1 microexon may be impacting neurogenesis by modulating cytoskeletal remodeling and associated signaling pathways. CLASP proteins, through stabilizing peripheral microtubules with the extracellular matrix, facilitate focal adhesion (FA) turnover by severing cell-matrix connections49,50. Given the connectivity of microexon-containing genes to the FAK-Src signaling pathway (Fig. 2g), we investigated whether the Clasp1 microexon modulates FA dynamics via FAK. FAK mediated signaling is achieved via initial phosphorylation of FAK at Y397, which promotes intermolecular interactions with Src-homology (SH)2 domain containing proteins51. This in turn results in the activation of downstream signaling events, including ERK phosphorylation. Phosphorylation of ERK subsequently promotes molecular changes associated with cell motility, differentiation or cell survival, depending on cellular context51,52. Interestingly, we observe an increase in FAK phosphorylation when the Clasp1 (A) microexon is completely skipped in ESCs and this increase becomes much more pronounced in differentiating neurons (Fig. 4g), resulting in increased ERK phosphorylation (Fig. 4h). In contrast, forced inclusion of the microexon has little effect on FAK phosphorylation in ESCs or neurons, although it does increase ERK phosphorylation in ESCs, and to a lesser extent in neurons (Fig. 4h), indicating a potential alternative signaling pathway activated by selective expression of the neuronal isoform. Taken together, these data reveal a previously unknown role for a Clasp1 microexon in modulating neurogenesis. In particular, the Clasp1(A) microexon-skipped isoform promotes FAK activation compared to the microexon-included isoform, which induces downstream ERK signaling, an event that is critical for the progression of neurogenesis53 (Fig. 4i). Similar to the loss of function of neurodevelopmental disorder genes such as SYNGAP154 and PCDH1255, we observe that perturbation of the Clasp1 (A) microexon results in the dysregulation of cytoskeletal signaling, as reflected through a significantly altered transcriptome and accelerated neuronal maturation.

Beyond the Clasp1 (A) microexon, we applied dCasRx-RBM25 to manipulate the splicing of microexons in Ptprf (A), Camta1, and Ptk2 (A) (Supplementary Figs. 6b and 6c). To investigate whether differential gene expression signatures arising from these perturbations correspond to neuronal differentiation gene expression programs, we correlated detected changes in DGE profiles with gene ranking by CellRank2. We observe significant positive correlations for Ptprf(A), Clasp1(A) and Gfra1 microexon exclusion and deletions, and negative correlations for Camta1 microexon exclusion, for both the bulk RNA-seq and CHyMErA-seq data (Supplementary Fig. 6d), confirming roles of these microexons in regulating neuronal gene expression timing. For the Ptk2 (A) microexon deletion, we observe a weak positive correlation (R=0.05, p=0.02; Pearson correlation) from the CHyMErA-seq data, possibly due to sparse data, and a significant negative correlation (R=-0.21, p=1.82×10−124; Pearson correlation) from bulk RNA-seq data. This suggests that this microexon may serve to promote neuronal differentiation, as part of an overall ‘push-and-pull’ mechanism involving different microexons within the same network of proteins that, collectively, function to promote or attenuate timing of different stages of neurogenesis. Overall, these data confirm the dysregulation of gene programs associated with neuron morphology and development observed in the CHyMErA-seq screen and extend these observations by highlighting disruption of downstream signaling cascades as a result of microexon perturbation.

Neuronal microexons control developmental timing of gene expression signatures linked to autism
We next asked whether the transcriptomic signatures detected in the CHyMErA-seq screen data overlap with those in the brains of individuals diagnosed with ASD. To address this, we utilized ‘Single-Cell rEgulatory Network Inference and Clustering’ (SCENIC)56 to define and score ‘regulon’ activity (i.e., data-derived co-expression networks of genes attributed to co-expressed transcription factors with proximal binding motifs) for individual cells. We focused our analysis on regulons that in a recent study57 show dysregulation in excitatory neurons in ASD compared to non-ASD brains. We compared scores for these ASD-dysregulated regulons between microexon-deleted cells and non-targeting cells, and observed significant differences in activation of Sox6, Ets2, Egr1, Elk1, Foxp1 and Sin3a regulons (P<0.05, Benjamini-Hochberg corrected two-tailed unpaired t-test) (Fig. 5a). Of these six dysregulated regulons associated with microexon deletion, three are downstream of MAPK-ERK signaling pathway (Ets2, Egr1 and Elk1), which has been implicated in core processes during neurogenesis58. Notably, in addition to the role of the Clasp1(A) microexon in augmenting MAPK-ERK signaling shown in this study, roles in MAPK signaling have been previously described for Bin159, Gfra143, Med2360–62, Ptprf63,64 and Ptk251 proteins (which also form a protein-protein interaction network, Fig. 2f) suggesting that the microexons in these proteins may convergently regulate this signaling pathway. Overall, we observe that 6 of 8 microexon deletion signatures show dysregulation of the Egr1 regulon, which is known to mediate the induction of p35, a neuron-specific activator of cyclin-dependent kinase 5 (Cdk5) downstream of ERK, an event that is essential for neurite outgrowth and neuronal differentiation65, thus highlighting a potential common pathway by which microexon-containing regulators of cytoskeletal processes regulate neurogenesis.

To further investigate whether genes implicated in ASD exhibit dysregulation in response to microexon deletion, we examined the transcriptional profiles of microexon deletions with premature differentiation phenotypes (i.e., Bin1, Clasp1(A), Gfra1, Ptprf(A), Med23, Ralgapb). A ‘module score’ representing expression of genes with genetic links to ASD (using ASD association evidence scores 1-3 from Simons Foundation Autism Research Initiative (SFARI)66) was derived for each single cell, and pair-wise comparisons of these scores were performed for microexon deleted cells and non-targeting cells in the neuroepithelial cell subcluster to investigate potential early activation of ASD-linked gene signatures. We observe a significant increase in the module scores for all deletions examined (P<0.001, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) indicating that skipping of these microexons promotes increased expression of these genes at an early developmental stage. To determine how these module scores deviate from non-targeted cells during development, we binned cells according to their calculated pseudotime and compared ASD-gene expression profiles between microexon-targeted (Bin1, Clasp1(A), Gfra1, Ptprf(A), Med23, Ralgapb microexons) and non-targeted cells (Fig. 5b). We observe that over pseudotime these premature maturation signatures include an early increase in expression of ASD-linked genes (pseudotime bins 3-6. p<0.05, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) followed by a lagging increase in non-targeting cells (pseudotime bin 7, p<0.05, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test), before reaching a plateau (Fig. 5c). This suggests that these microexons may tune expression of ASD-risk genes within neuronal transcriptomes during transient periods of neurogenesis. We also observe enrichment of ASD risk genes identified in three studies66–68 among the differentially expressed genes in response to perturbation of Clasp1(A), Ptprf(A), Ptk2(A), Camta1 and Gfra1 microexons in our bulk RNA-seq (Fig. 5d).

The convergence of the transcriptomic signatures derived from microexon deletion during in vitro neurogenesis on regulons and that are disrupted in the brains of individuals with ASD, and the altered timing of ASD-linked gene expression, further suggests a more widespread role of neural microexon dysregulation in neurodevelopmental disorders. In this regard, from an analysis of BrainGVEX data from the PsychENCODE consortium69 we note that several of the microexons linked to transcriptomic perturbations in our study show dysregulated splicing in the brains of individuals with Bipolar Disorder (Bin1, Ptk2(A), Ptk2(C), Ptprf(A)) and Schizophrenia (Bin1, Ptk2(A), Ptk2(C), Ralgapb) (P<0.05, Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test) (Supplementary Fig. 7). Taken together with the observations described above indicating functions for these microexons in controlling neuronal development, our results highlight a potential role for their dysregulation in contributing to the etiology of neurodevelopmental disorders.

Methods
Cell culture
CGR8 mouse embryonic stem cells (mESCs) were cultured on gelatin-coated plates in GMEM supplemented with 100 μM β-mercaptoethanol, 0.1 mM nonessential amino acids, 2 mM sodium pyruvate, 5000 units ml–1 of penicillin/streptomycin, 1000 units ml–1 of recombinant mouse LIF (all Life Technologies) and 10% ES fetal calf serum (Wisent). Cells were passaged using Trypsin (Life Technologies) and maintained at subconfluent conditions at 37 °C and 5% CO2. Cells were regularly monitored for mycoplasma infection.

In vitro neuron differentiation
In vitro differentiation of CGR8 mESCs into neurons was performed as previously described87. At day in vitro (DIV)-8, 5.1×106 cells were seeded on 15 cm non-tissue-culture-treated plates containing 30 mL of serum-free differentiation media (GMEM containing 5% KnockOut Serum replacement (Thermo Fisher Scientific), 0.1 mM nonessential amino acids, 1 mM sodium pyruvate, 5000 units ml–1 of penicillin/streptomycin and 100 μM β-mercaptoethanol). Cells were incubated for 3 days undisturbed. 50% media changes were conducted every 48 h. At DIV-5 the media was supplemented with 6 μM retinoic acid (RA) with subsequent media changes also containing 6 μM RA. At DIV0, cellular aggregates were dissociated with TrypLE Express (Gibco) for 5 minutes at 37 °C. Trypsinization was halted with soybean trypsin inhibitor (Thermo Fisher Scientific) and aggregates were homogenized by trituration with a 1 mL micropipette. Cells were pelleted for 5 minutes at 300xg. Neural progenitors were washed in N2 medium (Neurobasal-A medium (Gibco) with 1x N2 supplement: 2 mM glutamine and 5000 units ml–1 of penicillin/streptomycin (Gibco)), counted using a hemocytometer, and plated at a density of 1.5×105 cells/cm2 on poly-D-Lysine (Sigma-Adlrich) and laminin (mouse, Sigma-Aldrich) coated dishes. The plated neural progenitors were washed with N2 medium after 24 h to remove residual serum and non-adherent cells. At DIV2, the N2 medium was replaced with B27 medium (Neurobasal-A supplemented with antibiotics, 2 mM glutamine and 1xB27 supplement (Gibco)). Subsequently, differentiating neurons were provided with full medium changes with B27 on DIV4, and every 3 days onward.

Lentiviral transduction
HEK293T cells were seeded in DMEM (Wisent) supplemented with 10% FBS (Gibco) and 500 U mL-1 penicillin/streptomycin at a density of 5×105 cells per 6 well plate. 24 h later the cells were transfected with packaging vectors psPAX2 (Addgene: 12260), pMD2.G (Addgene: 12259) and the target vector at a 1:1:1 molar ratio. The next day the media was changed to high-BSA harvest media (DMEM supplemented with 1.1% (w/v) BSA and 5000 U mL-1 penicillin/streptomycin). 24 h later, virus containing media was collected, spun down at 500 rpm to pellet any packaging cells, and either immediately used for transduction, or snap-frozen in liquid nitrogen and stored at -80 °C. To generate the hgRNA library, transfections were scaled up to a 10 cm plate. Transductions were carried out by adding virus containing media and 8 μg mL-1 polybrene to the intended cells.

Cas12a variant vectors
Neomycin resistance was subcloned from the plenti-Lb-Cas12a-2xNLS vector (Addgene: 155046) into pRDA_112 (Addgene: 136475) via restriction cloning to allow for multiple drug selection of cells transduced with Lenti-Cas9-2A-Blast (Addgene: 73310). Addition of c-myc nuclear localization signals (NLSs) to the c-terminus, or removal of the n-terminal SV40 NLS of each vector, was achieved by site directed mutagenesis (Q5 Site-Directed Mutagenesis Kit (NEB)).

Cell line generation
CHyMErA cell lines. CGR8 cells were transduced with lentivirus carrying the Cas9-2A-BlasticidinR-expressing cassette (Addgene, no. 73310) and selected with blasticidin (6 µg ml-1) for 10 d. Cas9-expressing cell lines were then transduced with lentivirus carrying enAs- or Lb-Cas12a-2A-NeoR expression cassettes and selected with G418 (500 µg ml–1) for 14 d.

Gfra1 microexon-deletion clones. Cas9 guides targeting regions flanking the Gfra1 microexon or intergenic regions were cloned into px459 (Addgene: 62988) and transfected into CGR8s. Transfected cells were selected with 1 μg mL-1 puromycin selection for 3 days and plated by limiting dilution into 96-well plates for single-cell isolation. Isolated clonal cells lines were genotyped by genomic PCR of the edited region and analyzed by RT-PCR to confirm expression of the microexon excluded isoform.

dCasRx-RBM25 lines. Guide RNAs were designed either overlapping the microexon or within the downstream intron, to promote skipping or inclusion of the exon, respectively, and were scored using the TIGER algorithm88. Guide sequences were cloned into pLentiRNAGuide (Addgene: 138150) and transduced into CGR8 lines stably expressing dCasRx-RBM2511 (Addgene: 221002) and selected with 1 μg mL-1 puromycin for 3 days.

CROP-seq vector construction
Single Guides. The tracrRNA region of the original CROP seq-Guide-Puro (Addgene: 86708) was removed by site directed mutagenesis (Q5 Site-Directed Mutagenesis Kit). Guides were then cloned into this modified vector as previously described7 with minor modifications. Briefly, complementary oligonucleotides containing Cas9 and Cas12a spacer sequences separated by a buffer sequence were ordered, annealed and ligated into Esp3I-digested modified CROP-seq vector. To introduce the intervening Cas9 tracrRNA and Cas12a direct repeat sequences, a golden gate assembly was performed with Esp3I and a TOPO vector carrying Cas9 tracrRNA and either the As- or Lb-Cas12a direct repeat over 12 cycles ((1) 37 °C, 30 min; (2) 16 °C, 30 min; (3) 24 °C, 60 min; (4) 37 °C, 15 min; (5) 65 °C, 10 min; steps 1–3 were repeated for 11 cycles) using a vector/insert ratio of 1:25. Ligations were transformed into Stbl3 chemically competent cells (Invitrogen) and successful transformants were confirmed by Sanger sequencing.

Screen library. For construction of hgRNA libraries, Cas9 and Cas12a gRNA sequences were cloned into a lentiviral vector via two rounds of cloning (as in ref. 7) with slight modifications. Oligo pools (Twist Biosciences) were designed carrying 20-nt Cas9 and 23-nt Cas12a guide sequences separated by a 32-nt buffer sequence flanked by Esp3I restriction sites, all flanked by short sequences containing BfuAI restriction sites. Oligonucleotides were amplified by PCR over ten cycles using NEB Q5 Ultra polymerase (NEB) ((1) 98 °C, 30 s; (2) 98 °C, 10 s; (3) 53 °C, 30 s; (4) 72 °C, 10 s; (5) 72 °C, 2 min; steps 2–4 repeated for nine cycles). Amplified oligos were purified on a PCR purification column (GeneJET PCR purification kit; Thermo Fisher Scientific), and an aliquot was run on a 2% agarose gel to check purity. The modified CROP-seq vector backbone was digested with Esp3I (NEB) overnight at 37 °C. The digested backbone was dephosphorylated with rSAP (NEB) for 1 h at 37 °C and gel purified (GeneJET gel extraction kit; Thermo Fisher Scientific). Amplified oligos were digested with BfuAI (NEB) and ligated into the digested modified CROP-seq backbone using T4 ligase (NEB) in a combined reaction overnight at 16 °C using a 3:1 insert to vector ratio. The purified ligation reaction was used to transform Endura competent cells (Lucigen) by electroporation (1-mm cuvette, 25 uF, 200 Ω, 1600 V), and a sufficient number of cells were plated on 15-cm ampicillin Luria–Bertani (LB) agar plates to reach a library coverage > 1000-fold. Bacterial colonies were scraped from the plates, pooled and pellets were collected. The Ligation 1 library plasmid was extracted and purified (GeneJET Maxi-Prep Kit; Thermo Fisher Scientific).

In the second step, Cas9 tracrRNA and the Cas12a direct repeat were inserted into the pooled library. The Ligation 1 plasmid library was digested overnight using Esp3I (NEB) and dephosphorylated using rSAP (1 h, 37 °C) and purified on a PCR purification column. A TOPO vector carrying Cas9 tracrRNA and the Cas12a direct repeat was digested using Esp3I and subsequently ligated into the digested pLCHKO-Ligation 1 vector overnight over 12 cycles ((1) 37 °C, 30 min; (2) 16 °C, 30 min; (3) 24 °C, 60 min; (4) 37 °C, 15 min; (5) 65 °C, 10 min; steps 1–3 were repeated for 11 cycles) using a vector/insert ratio of 1:25. The purified ligation reaction was used to transform Endura competent cells by electroporation (1-mm cuvette, 25 uF, 200 Ω, 1600 V), and a sufficient number of cells were plated on 15-cm ampicillin LB agar plates to reach a library coverage of > 1000-fold. Bacterial colonies were scraped from the plates, pooled and pellets were collected. The Ligation 2 library plasmid was extracted and purified (GeneJET Maxi-Prep Kit; Thermo Fisher Scientific).

CHyMErA efficiency testing
CGR8 mESCs were transduced with lentivirus carrying an hgRNA-expressing modified CROP-seq vector targeting microexons or a non-targeting control. The day after transduction, cells were subjected to 1 μg ml−1 puromycin selection for 3 days to ensure complete selection. Genomic DNA (gDNA) from cells was purified (PureLink Genomic DNA Mini Kit; Themo Fisher Scientific) and the targeted region was amplified by PCR. PCR products were resolved on agarose gels and band intensity was quantified with ImageJ. Percent editing was calculated as the quotient of the size normalized band intensities:
% Deletion = 100 * (Deletion Band intensity / Deletion Band length) / (Uncut Band intensity / Uncut Band length + Deletion Band intensity / Deletion Band length) (1)

To independently determine approximate editing efficiencies of the eCHyMErA method, amplified PCR products shown in Fig. 1c were analyzed by Oxford Nanopore long-read sequencing (performed by Plasmidsaurus Inc.). Raw reads were aligned to sequences corresponding to the amplified regions. Editing efficiency was quantified by comparing the relative proportions of reads mapping to Cas9/Cas12a deleted regions versus adjacent unedited sites. Editing efficiencies estimated by quantification of gel band intensities and long-read sequencing were within 10% of each other.

Library validation
gDNA from transduced CGR8s was harvested and integrated hgRNA sequences were amplified using Q5 Ultra polymerase (NEB) ((1) 98 °C, 30 s; (2) 98 °C, 10 s; (3) 65 °C, 30 s; (4) 72 °C 20 s; (5) 72 °C 2 min; steps 2-4 were repeated 25 times) primers:
CROP-Cas9-F_s[N]: ACACGACGCTCTTCCGATCT[N]CTTGTGGAAAGGACGAAACACC
CROP-Cas12a-R_s[N]: CTGGAGTTCAGACGTGTGCTCTTCCGATCT[N]GTGTCTCAAGATCTAGTTACGCCAAGC
Where [N] indicates intervening nucleotides to generate staggered amplicons to increase library diversity.

These amplicons were then subject to a second round of PCR ((1) 98 °C, 30 s; (2) 98 °C, 10 s; (3) 55 °C, 30 s; (4) 72 °C 20 s; (5) 72 °C 2 min; steps 2-4 were repeated 14 times) to add sequencing adapters. The resulting amplicons were gel purified and paired-end sequenced on an Illumina MiSeq instrument using the following cycling parameters: 22 dark cycles, read 1: 27 cycles, index 1: 8 cycles, index 2: 8 cycles, 35 dark cycles, read 2: 27 cycles. Raw fastq files were aligned to an index of hgRNA sequences using Bowtie89 to determine library representation.

Single nuclear RNA-seq library preparation
DIV7 neurons were harvested with Hank’s Balanced Salt Solution (Gibco) supplemented with 20 Uml−1 Papain (Worthington), 5 mM MgCl and 5 μg mL−1 DNase I (Worthington), counted and pelleted and immediately snap-frozen in liquid nitrogen and then stored at −80 °C. Nuclei extraction and fixation were performed as previously described26, except for the use of a modified CST lysis buffer90 plus 1% of SUPERase In RNase Inhibitor (AM2696). Nuclei quality was checked with DAPI and Wheat Germ Agglutinin (WGA) staining. Sci-RNA-Seq3 libraries were generated as previously described91 using three-level combinatorial indexing. The final libraries were sequenced on an Illumina NovaSeq with read cycles as follows: read 1: 34 bp, read 2: 69 bp, index 1: 10 bp, index 2: 10 bp.

Hemi-nested hgRNA enrichment PCR
Enrichment PCR was adapted from92 for capture of Cas9 and Cas12a spacer sequences. hgRNA sequences were sub-amplified from the barcoded library using three rounds of PCR. In the first round, primers annealing upstream of the hgRNA sequence and to the P5 barcode were used to amplify the region of the guide-barcode transcript. In the second round, primers annealing upstream of the Cas9 and Cas12a sequences were added to sub-amplify the spacer sequences and provide handles for Illumina adaptors. Illumina adaptors were added in a final round of PCR. After each PCR reaction, libraries were purified with Ampure XP beads (Beckman Colter), and product sizes were confirmed by agarose gel electrophoresis. Final libraries were sequenced on an Illumina NovaSeq with the read cycles: read 1: 34 bp, read 2: 69 bp, index 1: 10 bp, index 2: 10 bp.

Immunofluorescence and neurite length analysis
Gfra1 microexon deletion clones or wild-type CGR8 mESCs were differentiated into neurons (see In vitro Neuron Differentiation). At DIV0 neural progenitors were plated on Laminin (Sigma-Aldrich) and Poly-D-Lysine (Sigma-Aldrich) coated coverslips in 12-well plates. At DIV3 and DIV7, cells were washed with PBS and fixed in 4% PFA for 15 min. Cells were washed with PBS and permeabilized in 0.25% Triton X-100 for 10 min. Cells were then washed 3 times in PBS, blocked in PBS with 2.5% BSA for 1 h, then incubated overnight with 1:2000 chicken anti-MAP2 antibody (Abcam; 5392) in PBS with 2.5% BSA. The next day cells were washed 3 times with PBS and incubated with goat anti-chicken IgG Alexa-Fluor 647 conjugated secondary antibodies and 1:1000 DAPI in PBS with 2.5% BSA for 1.5 h. Slips were washed and mounted on glass slides and imaged on an Axiovert 200m inverted fluorescence microscope (Zeiss). Image IDs were randomized to blind the researcher to the cell type analysed, and neurites from MAP2 stained neurons were quantified with the NeuronJ plugin93 for ImageJ.

Immunoblotting
For analysis of FAK and ERK phosphorylation, CGR8 mESCs or DIV7 neurons were washed in ice-cold PBS and lysed in ice-cold RIPA buffer (10 mM Tris-HCl, 150 mM NaCl, 0.1% SDS, 0.5% sodium deoxycholate, 1% Triton X-100, 1 mM EDTA) supplemented with cOmplete protease inhibitor cocktail (Roche) and PhosSTOP phosphatase inhibitor (Roche). Lysates were collected in 1.5 mL Eppendorf tubes, sonicated, and protein content was determined using the Bradford assay (BioRad). Samples were mixed with Laemmli buffer and 25 μg of protein was resolved on a polyacrylamide SDS-PAGE gel before overnight transfer to a PDVF membrane. Membranes were blocked for 1 h in either 5% milk or 5% BSA in TBS-T for endogenous and phospho-antibodies, respectively. Membranes were then washed 3 times in TBS-T and incubated with 1:1000 anti-FAK (Cell Signaling, 3285) or 1:1000 anti-p44/42 MAPK (Erk1/2) (Cell Signaling, 9102) in 5% milk TBS-T, or 1:1000 anti-Phospho-FAK (Y397) (Cell Signaling, 3283) or 1:1000 anti-Phospho-p44/42 MAPK (Erk1/2) (Thr202/Tyr204) (Cell Signaling, 9101) in 5% BSA TBS-T overnight at 4 °C. Membranes were then blotted for tubulin (Sigma-Aldrich, T-6074) as a recovery and loading control. Detection was achieved using anti-rabbit/mouse IgG HRP-linked secondary antibodies (Invitrogen) and chemiluminescence (PerkinElmer, NEL120001EA).

RT-PCR and RNA-seq
Total RNA was extracted using the RNeasy Mini Kit (QIAGEN) as recommended by the manufacturer. For analysis of microexon splicing, RT-PCR was performed with the OneStep RT-PCR kit (Qiagen) and products were resolved on 3.5-4% agarose gels stained with ethidium bromide. Exon inclusion and exclusion were quantified by ImageJ. For RNA-seq, RNA was similarly isolated and treated with on-column DNAse (Qiagen). Libraries were then generated for polyA+ stranded RNA-seq with the NEBNext Ultra II Directional kit and 150 nucleotide paired-end reads were sequenced on an Illumina NovaSeq 6000 at an average depth of 25 million reads.

Computational procedures
Microexon-targeting library design. For the microexon targeting hgRNA library, mouse microexons were selected according to the following criteria. 1) they correspond to an orthologous human exon. 2) they exhibit strong differential inclusion during in vitro neuronal differentiation30, with ≥50% increase in inclusion level between undifferentiated cells (DIV-8) and neurons (DIV7-28), (3) they are spliced in transcripts with expression ≥5 corrected Reads Per Kilobase-pair and Million mapped reads (cRPKM), and (4) they exhibit differential inclusion in an in vivo analysis of mouse brain development31. Microexons from this set were further selected on the basis of their location in genes with defined roles in neurogenesis, and that have been implicated in ASD, using ASD association evidence scores 1-3 from the Simons Foundation Autism Research Initiative (SFARI).

From this set of microexons we extracted and scored Cas9 (CRISPOR32) and Cas12a (CHyMErA-NET7) guides in the flanking intronic regions, ensuring minimally a 150bp distance from annotated upstream and downstream exon splice sites, and a 100bp distance from the microexon splice sites. The highest scoring pairs of guides were used to direct hgRNA library design. These selection steps resulted in 37 microexons with introns that could accommodate multiple pairs of upstream and downstream targeting guides. Intronic guides were designed in a similar fashion, except that both Cas9 and Cas12a guides were restricted to either the upstream or downstream intron, leaving a distance of at least 300bp from any splice site to ensure that the splicing of the microexon or neighboring exons would not be affected. Guides targeting constitutive exons were designed using CRISPICK94,95, and non-targeting and intergenic control guides were derived from a previously published library7.

Single nuclear RNA-seq data analysis
Raw data processing. Raw sequencing reads were demultiplexed based on i5/i7 PCR barcodes and processed using a modified pipeline from26 implemented in a Snakemake framework. First, barcodes and unique molecular identifiers (UMIs) were extracted from read 1 of FASTQ files followed by trimming of read 2. Reads were then aligned to the mouse genome (mm10) with STAR short-read aligner (v2.6.1d) using Gencode vM12 gene annotations and quality filtered. Reads were deduplicated based on UMIs, and summarized into a count matrix of M genes x N nuclei, which was converted into AnnData object in Scanpy (v1.9.3)96.

Quality filtering and normalization. Doublets were removed with Scrublet97 and cells with a detected number of genes, total RNA counts, or mapping rate greater than five absolute deviations from the median were discarded, as were cells with > 1% mitochondrial RNA. Total counts per cell were normalized and log transformed.

Dimensionality reduction and clustering. The top 2000 most variable genes were extracted and used in principal component analysis (PCA) to reduce the dimensionality of the data. Neighborhood graphs were constructed using 30 principal components and data were visualized with Uniform Manifold Approximation and Projection (UMAP). Clustering was performed using the Leiden algorithm (leidenalg v0.10.2) with a resolution of 0.1 to identify broad cell clusters and 0.8 to identify subclusters. A specific subcluster (‘15’) was noted to have aberrantly high RNA counts, and expression signatures consistent with doublets, and therefore was excluded from downstream analysis.

Clusters were manually annotated based on expression of known cell type marker genes and label transfer using an embryonic mouse dataset and named accordingly (see Cell-type label transfer).

Diffusion map and pseudotime analysis. A diffusion map was computed using sc.tl.diffmap, and a root cell was identified based on the minimum value in the third diffusion component corresponding to cells with transcriptomes indicative of early stages of differentiation. Pseudotime analysis was performed using the diffusion pseudotime (DPT) method (sc.tl.dpt)98.

Cell-type label transfer. Preprocessed data were loaded into a Seurat (v4.3.0)99 object. Reference data from Cao et al. (ref. 26) or Pijuan-Sala et al. (ref. 100) underwent similar preprocessing, including normalization, scaling, and dimensionality reduction via PCA. To transfer labels from reference data, integration of the query and reference datasets was performed using Seurat’s anchor-based label transfer approach. Anchors were identified between the datasets using the first 30 principal components, and cell-type annotations from the reference dataset were transferred to the query dataset.

CellRank2 analysis. To explore neurodevelopmental trajectories and cellular dynamics, pseudotime analysis coupled to CellRank2 was employed to infer the progression of cells along a continuous developmental path. Transition matrices were computed using the PseudotimeKernel from CellRank2 (v2.0.4)33,42, followed by embedding projections. The macrostate analysis and terminal state probabilities were calculated with the ‘Generalized Perron Cluster Cluster Analysis’ (GPCCA) estimator (pygpcca v1.0.4)101, and coarse-grained transitions were visualized. Lineage driver genes were identified and analyzed by computing their correlation with specific lineage states. Upregulation of neuronal lineage ‘driver genes’ was calculated by two-sided Wilcoxon rank sum test comparing up- or down-regulated gene ranking (FDR<0.05) with background genes that were not differentially expressed, for either single nuclei or bulk RNA-seq.

Guide assignment. Sequencing data from hgRNA libraries were processed using a Snakemake pipeline. Raw FASTQ files were processed to attach UMIs using custom scripts adapted from ref. 26, incorporating ligation and reverse transcription barcodes. Reads were aligned to a non-redundant hgRNA index using Bowtie102. Aligned reads were filtered to retain high-quality mappings. hgRNA-UMI pairings were considered chimeric if they represented less than 80% of the reads for a given hgRNA-UMI pairing within a given cell, and were removed using a custom Python script, and the resulting non-chimeric BAM files were indexed. Deduplication was performed using UMI-tools103 to retain unique molecules.

Aligned reads with UMIs were processed to assign guides to individual cells. Reads were grouped by cell and guide, and UMIs were counted to determine the number of unique reads per cell and guide. Guides were assigned to cells by prioritizing the guide pair with the highest combined UMI counts or, for non-redundant guides, where possible, the single guide with the highest UMI count. Perturbation targets were annotated based on guide metadata, and final assignments were compiled.

Perturbation calling
Normalization to non-targeting cells. To normalize individual perturbation signatures to cells at a similar stage in development, gene expression data was normalized to the average expression in the k=30 nearest non-targeting cells (using a similar strategy described in ref. 35). Principal component analysis (PCA) was performed, and the 30 nearest neighbor non-targeting cells were calculated for each cell based on Euclidean distance from the top 40 PCs. We calculated a normalized gene expression vector for each cell xi subtracting the average expression of the neighboring cells Ni for each gene:
xei = xi - xNi (2)
where:
xNi = 1/k * sum(xj) (3)

Perturbation assignment. Raw count data for each perturbation were aggregated into pseudo-bulk profiles by cluster and guide RNA assignment. Differentially expressed genes were identified by pairwise comparison of pseudobulk expression profiles for cells expressing guides targeting the same microexon, and non-targeting cells. For each comparison, we used the likelihood ratio test (LRT) from edgeR (v3.38.4) to calculate log (fold change) and FDR-adjusted p values. Only pseudobulks comprising at least 5 aggregated cell expression profiles, and perturbations with at least 25 total cells were considered.

For each perturbation, the normalized expression values calculated above (Normalization to non-targeting cells) corresponding to the DEGs (minimum of 5 differentially expressed genes; FDR<0.05) were used as input for classification.

To classify cells as perturbed or non-targeting we applied quadratic discriminant analysis (QDA), implemented in sci-kit learn (v1.2.2), to calculate classification assignments based multivariate Gaussians fit on normalized expression values xei. QDA models each class as a multivariate Gaussian with its own class-specific covariance matrix Σk and mean μk and assigns classifications using Bayes’ rule:
P(yi=k | xi) = P(xi | yi=k) * P(yi=k) / P(xi) (4)
for which we calculate the class conditional likelihood as:
P(xi | yi=k) = (2π)^d/2 * |Σk|^-1/2 * exp(-1/2 * (xi-μk)^T * Σk^-1 * (xi-μk)) (5)

Priors were inferred from class proportions in the data. QDA was performed across 1000 bootstrapped replicates. During each bootstrap iteration, cells were randomly resampled with replacement from the perturbed and non-targeting groups to fit the model and predictions were made for all cells. The bootstrapped classifications were aggregated to compute an overall classification accuracy for each cell. Statistical significance of these classifications was assessed using a permutation-based approach. Labels for perturbations were shuffled to create a null distribution of classification accuracies. For each cell, an empirical p value was calculated by comparing its observed classification accuracy to the null distribution:
p = (1 + sum(I(p^i <= pnull,n))) / (N + 1) (6)

where p^i is the observed classification accuracy of cell i, and pnull,n is the accuracy in the nth permutation.

These p values were corrected for multiple comparisons (Benjamini-Hochberg), and cells with a corrected p-value < 0.05 were classified as ‘knockout’, and cells with p > 0.05, or perturbations with fewer than 25 cells were classified as ‘non-perturbed’.

Guide correlation analysis. To investigate how transcriptional signatures between guides correlate compared to non-targeting cell populations, we compared pseudobulks of DEGs (FDR<0.05) for each microexon deletion. We computed a Pearson correlation between log-normalized expression values from cells expressing guides targeting the same microexon, and those expressing non-targeting, or intergenic targeting guides. Significance was assessed by Benjamini-Hochberg corrected two-sided Wilcoxon rank sum test.

Comparison of microexon-deletion and gene knockout cells. For each microexon deletion that resulted in a detectable transcriptomic phenotype, we identified differentially expressed genes (FDR<0.1) and performed principal component analysis (PCA) on their normalized expression values across single cells (see ‘Normalization to non-targeting cells’ in the Methods). Cells were embedded in PCA space using the top 40 principal components, or fewer if limited by sample size. To assess transcriptomic differences, we performed pairwise comparisons using Hotelling’s T² test between microexon-targeting cells and those expressing microexon deletion, gene knockout guides, or intergenic negative control guides, relative to non-targeting controls. This approach enabled statistical testing of whether transcriptomic signatures induced by microexon deletion were distinct from those induced by gene knockout or negative control perturbations.

Temporal analysis of differential gene expression. Gene signature scores were computed per cell by summing normalized log-expression values (log1p-transformed) of genes significantly upregulated in Gfra1 microexon deletion (KO) versus non-targeting (NT) cells (logFC>0, FDR<0.05). Scores were z-normalized across all cells. A generalized additive model (GAM) was then fitted to the scaled scores as a function of pseudotime using a smooth interaction with guide type. Model fits and 95% confidence intervals were computed using the mgcv (v1.9_3) package in R104.

Gene ontology enrichment and overlap analysis
GO enrichment was performed using g:Profiler (v0.2.3) implemented in R105. Differentially expressed genes were used as input with a custom background set of genes comprising all genes involved in DEG testing for each pair-wise comparison.

To assess concordance between GO enrichments for CHyMErA-seq and bulk RNA-seq (adjusted P < 0.05), we used a contingency table of term overlap. Fisher’s exact test was applied to determine the statistical significance of overlap and odds ratios from both platforms were used to compute combined enrichment scores. GO terms were reduced using semantic similarity matrices using the rrvgo106 (v1.18.0) package in R, with grouping based on parent terms at a threshold of 0.8.

Gene set enrichment analysis (GSEA)
GSEA was performed using the fgsea package in R107. For analysing single cell expression signatures, genes were ranked by their log2(fold-change) values from edgeR. Gene sets were retrieved from MSigDB108,109, filtered for GO terms comprising between 15 and 150 genes. For analysing the CellRank driver genes, genes with a significant value (p<0.05) were ranked according to their correlation with the neuronal trajectory, and gene sets were filtered for GO terms comprising between 15 and 500 genes.

SCENIC analysis
Gene regulatory network (GRN) inference and downstream analysis were performed using the pySCENIC pipeline (v0.12.1)56,110 using preprocessed data (see Single Nuclei RNA-seq Data Analysis). For regulatory network analysis, the gene expression data were converted into a Loom file format, which was then used as input for GRNBoost2111, to infer gene regulatory networks. After generating the adjacency matrix of inferred gene regulatory interactions, regulon pruning was performed to reduce redundancy and focus on co-expression of genes with transcription factors with proximal binding motifs. AUCell110 was then used to calculate regulon activity for each single cell.

Regulons previously identified as being dysregulated in the brains of individuals with ASD57 were analysed for differences in average scores between microexon perturbed cells and non-targeting control cells. Regulon scores were z-normalized and p values were calculated for each pairwise comparison by two-tailed unpaired t-test followed by Benjamini-Hochberg correction.

ASD gene-expression timing analysis
Expression modules for individual cells in the neuroepithelial cluster were calculated using Seurat’s AddModuleScore function with ASD-linked genes in the SFARI gene dataset that were differentially expressed in at least 1 microexon deletion (FDR<0.05) as the input genes. These scores were then compared with a population of ‘non-targeting’ cells in the same cluster and significance was assessed by Wilcoxon rank-sum test. To test the developmental timing of ASD-linked gene dysregulation, cells from microexon deletions with premature differentiation phenotypes (Bin1, Clasp1(A), Gfra1, Ptprf(A), Med23, Ralgapb microexon deletions) or non-targeting cells, were binned according to pseudotime and pairwise comparisons of SFARI gene module score for each pseudotime bin were calculated by Wilcoxon rank-sum test. P values were adjusted by Benjamini-Hochberg correction.

ASD gene enrichment analysis
Enrichment was calculated for genes from ASD datasets in up- and downregulated genes (FDR<0.05) for each microexon exclusion directed by dCasRx-RBM25, or for clonal line deletion of the Gfra1 microexon. Odds ratios were calculated by Fisher’s exact test and p values were adjusted by Benjamini-Hochberg correction.

Bulk RNA-seq differential gene expression analysis
Clonal cell lines harboring Gfra1 microexon deletions, or simultaneously prepared clonal cell lines without genetic perturbation, were differentiated into neurons in vitro (as above), RNA was extracted at DIV7, and RNA-seq was performed (see RT-PCR and RNA-seq methods above). For dCasRx-RBM25 lines, cell lines expressing microexon inclusion or exclusion gRNAs were differentiated into neurons in vitro and RNA from two biological replicates for mESCs or DIV7 neurons was collected and RNA-seq performed. Gene-level quantification was performed with Salmon (v1.10.3)112. Differentially expressed genes were calculated with the quasi-likelihood F-test in edgeR (v3.38.4) comparing Gfra1ΔMIC/ΔMIC lines with wild type clones or samples expressing microexon exclusion guides compared to microexon inclusion guides for dCasRx-RBM25 lines.

PPI network analyses
We used the STRING database (v12.0) to identify potential physical interactions among 145 mouse proteins that contain neuronal microexons exhibiting switch-like inclusion both in vivo and in vitro during neurogenesis (see Microexon-targeting library design). To construct the interaction network, we included all detected protein-protein interactions with at least medium confidence (score ≥ 0.4). We then clustered the network using the Markov Cluster Algorithm (MCL) implemented in Cytoscape, applying a granularity parameter of 2 to balance cluster specificity and connectivity. We performed functional enrichment analysis on each cluster using STRING.

Microexon quantification from PsychENCODE
We quantified inclusion levels of microexons from RNA-seq collected from human brain samples available through PsychENCODE69. Microexon quantification of BrainGVEX bulk RNA-seq samples data was performed using MicroExonator31 with the following parameters: filter_mode = “unbiased”, use_corrected_PSI = True, min_PSI=0.1, and min_reads_PSI=5. Novel microexon discovery was skipped (skip_discovery=True), and quantification was limited to microexons annotated in VastDB3. Microexon inclusion levels were aggregated by primary diagnosis, and statistical comparisons were performed using the two-sided Wilcoxon rank sum test.

Data availability
Raw and processed data associated with the CHyMErA-seq screen have been deposited in GEO and can be found under the accession number (GSE291610). Raw bulk RNA-seq data and quantifications analysed in this study have been deposited in GEO and can be found under the accession number (GSE291608). Processed bulk RNA-seq data analysed in Fig. 1f can be accessed from VastDb (https://vastdb.crg.eu/3). Source data are provided with this paper.

Code availability
All code and processed data required to produce the main figures can be found on GitHub (https://github.com/stevendupas/chymera-seq) (https://doi.org/10.5281/zenodo.17692659).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Dupas and colleagues present CHyMErA-seq by extending the CHyMErA screening method developed earlier by the group. CHyMErA-seq allows one to delete neural microexons in a high-throughput manner using CRISPR/paired gRNAs and evaluate the functional consequence of the genetic perturbation in individual cells by split-pool single nuclei RNA-seq. This strategy is used to test a set of 37 conserved and neuron-specific microexons in mESCs during in vitro neurogenesis. The study finds neuronal transcriptomic shifts caused by microexon deletion and further examines the functional roles of two microexons (Gfra1 and Clasp1) in neurodevelopment. Deletion of microexons has also been implicated in altering expression of ASD-risk gene networks, providing insights into how dysregulated microexon splicing may contribute to disease. The functional evaluation of dynamic alternative splicing patterns, especially high-throughput technologies, is highly valuable. Indeed, within this class, the biological requirements of microexons are particularly important due to their potential connection with neurodevelopmental disorders.\nWhile the authors have made great strides in the execution of this high throughput strategy, this reviewer has some concerns regarding the approach itself, and the robustness of the downstream gene expression changes, and their relevance to maturation and disease. In particular, the manuscript will benefit from strengthening key aspects, including microexon deletions, as the entire body of work depends on this strategy, characterization of the approach, and phenotyping of neurons with exon deletions.\nCHyMErA-seq screening:\n1. Microexon deletion: Figure 1, the authors show %del to account for deletions. This is likely an overestimation as smaller PCR amplicon and typically favored. Even if one ignores this bias, the estimated deletion frequency (e.g., top panels in Figure 1d) does not seems to be consistent with the gel images in some cases.\nMoreover, the level of deletion in a bulk sample is not an especially useful metric, as what is more relevant is %cells with homozygous deletions. Carefully evaluating (for example with ddPCR), and reporting these numbers will help the audience understand the robustness of this dual Cas approach. The same problem applies to the application of CHyMErA-seq, where it is not clear to what levels homozygous microexon deletions occurs.\n2. The representation of guide pairs in the sequenced nuclei library (and potential sources of bias) needs further illustration. The authors have targeted 37 exons with 9 guide pairs each. How many cells were recovered per guide pair for each target that passed QC should be included. How is the correlation of DGEs between guide pairs targeting the same microexon affected by cell#?\n3. A major caveat of the CHyMErA-seq approach, when applied to mESCs, is the low editing efficiency, and cells with microexon deletion cannot be isolated directly. The use of QDA to classify cells as KO vs non-perturbed based on gene expression profile is a creative way to deal with the KO vs non-perturbed issue. However, this is an indirect workaround, and the reliability of the approach remains to be determined.\nThe authors compared CHyMErA-seq derived gene expression changes in response to microexon deletion and that determined by bulk RNA-seq of clonal samples, and the correlation is quite moderate (Figure 4d,e and Extended Data Figure 5d), which is very concerning.\n4. It seems that no replicate screening was performed, which adds additional challenge to evaluate the robustness of the screening.\n5. In some parts, the description of bioinformatics analysis is not very clear. For example, it is difficult to figure out how perturbation calling was done precisely. The authors might consider using less narratives and include some mathematical descriptions.\nPhenotyping:\n6. Are the 37 microexons frame preserving? This reviewer assumes so, but it will be helpful to clearly indicate this.\n7. The authors have positive control gRNA pairs that target constitutive exons, but do not provide much detail on the transcriptomic characterization of these cells.\n8. Pseudotime analyses indicates that the KO cells are broadly enriched in the neuroepithelia clusters (Fig 2d). This suggests there are differentiation defects and KO cells are stuck in the neuroepithelium state? As the microexons are typically not included to significant levels at this stage of differentiation, how do the authors explain KO cells enrichment in neuroepithelia. Moreover, this seems incongruent with the idea that microexons deletions increase maturation.\n9. More can be done to solidify/validate gene expression changes observed in the snRNA data. While the authors decided to focus on Gfra1 and Clasp1 for additional analyses and made efforts to generate clonal cell lines with microexon deletion (or skipping/inclusion using dCasRx-RBM25). The evidence provided to support the functional role of each exon is less than compelling in the current form. For example, for Gfra1, the transcritomic changes inferred from the screen essentially has no correlation with that from analysis of clonal cell lines (Extended Data Figure 5d). For Clasp1, the impact of microexon on downstream phosphorylation appears to be quite moderate.\n10. The authors used a KS-test to evaluate alterations of gene expression correlation with neuron differentiation rank. It is quite difficult to understand what these significant differences really mean.\n11. Ptk2 microexons are known to dramatically elevate autophosphorylation. Do the authors detect signatures of modulated Ras signalling/other downstream signaling activity? Similarly, do the authors observe elevated morphology or migration changes downstream of Ptk2 microexon modulation?\n12. GDNF signaling is known to promote survival and growth of neurons. The authors find that Gfra1 microexon deletion modestly increases neurite growth and stimulates expression of mature neuronal genes. Is this a consequence of increased GDNF signalling? i.e., does microexon inclusion decrease Gfra1 signaling? Can the authors supplement neurons with GDNF to amplify the neuritogenesis phenotypes?\n13. Camta1 microexon does not appear to be included in mESCs. Nevertheless, the changes in gene expression upon deletion in mESCs appears to enrich downregulated genes that are affected in autism (Fig 5d). this feature appears significant only in the mESCs and not in neurons. Similar results are also obtained for Clasp1 microexon as well. How do the authors explain observations of transcriptomic phenotypes in mESCs that do not express much of the microexon, but not in neurons from the same founders?\nOther minor points:\n14. Explain hgRNA\n15. Line 689, missing “)” after CRISPOR.\n16. Line 755, paring --> pairing\n17. Regulon scores were z-normalized and p values **and** were calculated for each pairwise comparison by two-tailed unpaired t-test followed by Benjamini-Hochberg correction.\n18. Fig. 2b, colors for some of the clusters are too light to be visible\n19. Fig. 2e, what are the columns?\n20. Fig. 4 legend, “Empirical cumulative distribution functon(ECDF) of differentially expressed genes (DEGs, FDR < 0.05) associated with microexon inclusion compared to **inclusion** for mESCs.”", "reviewer_2": "This manuscript presents CHyMErA-seq, a platform designed to profile transcriptomic changes at single-cell resolution in response to targeted splice isoform perturbations. Applying this approach to 36 neuronal microexons, the authors identify consequences of these alternative splicing events on gene expression dynamics during neurogenesis. Several microexon-containing genes tested here are implicated in neuron projection development and converge on signaling pathways such as FAK-Src. Functional analyses reveal that specific microexon deletions studied here (similar to the effect of microexons studied elsewhere) lead to increased neurite outgrowth and in general a premature activation of neuronal gene expression programs. The authors also observe parallels between microexon-dependent transcriptomic signatures and gene expression patterns observed in autism spectrum disorder, in line with previous studies.\nThis study builds on a large body of research from the Blencowe and other labs, which provided the basis for the experimental tools applied here. In the present study, the authors improved editing efficiency and coupled the CHyMErA pipeline with scRNA-seq readout. In addition, most if not all the main results obtained using scRNA-seq here confirm or are in line with previous studies, as the authors themselves point out repeatedly (activation of neuronal gene expression signatures, increased neurite outgrowth, link to ASD,…). On the one hand, this is an encouraging proof-of-concept of the experimental approach, but on the other hand raises questions about the novelty of the findings. In my opinion, it will be necessary to delineate more specifically which results are novel, which are confirmatory, and where the single-cell and isoform-level perspective (which is the main “selling point” of this study) allowed new biological insights into a previously studied biological question.\nSpecific points and suggested analyses are detailed below.\nEspecially considering the necessity of pointing out the novelty of the findings compared to previous studies, it would be advisable to be more specific in the wording of the results and avoid vague statements such as the ones below:\nPage 12: “Overall, these data confirm transcriptomic signatures detected in the CHyMErA-seq screen and further support important new functions of microexons in neurogenesis”\nPage 14: “Collectively, these findings highlight the importance of isoform-specific protein functions of genes in the control of critical and dynamic processes during neurogenesis, and their dysregulation in neurodevelopmental disorders.”\nPage 15: “Applying this system to neuronal microexons revealed unexpected functional properties [which ones?] of this important class of alternative splicing in the context of neurogenesis, as well as new insights into the multilayered links between microexons and brain disorders [I do not agree that the manuscript contains enough novel insights in this respect to justify this sentence].”\nSpecific comments:\nPage 8/Fig 2: “The Ptk2 microexons (denoted as ‘A’, ‘B’ and ‘C’) display coordinated splicing patterns in vivo…”: Their splicing is coordinated, but here the authors delete one of them at a time: Consequently, deletion of microexon alone will create an artificial protein that is not usually produced (because they are included or excluded together). Is there in vivo evidence for the existence of transcripts with two of the three microexons included and one skipped? In addition, the authors state that “notably, we observe transcriptomic signatures for all three microexons targeted in the Ptk2 gene (Fig. 2e)”. First of all, this should be rephrased to be less vague and more accurate: (a) …transcriptomic signatures of the DELETION of all three microexons…; (b) …we observe DISTINCT transcriptomic signatures…; (c) these apparent differences between the effect of the deletion of the three microexons should be discussed (in particular, B and C are located in the same domain, but yet deletion of each of them shows very different effects on the transcriptome).\nFig 3a and b: Figure legend and description in the text is very minimalistic. It is not clear to me if the focus here is on gene expression changes between Gfra1 microexon knockout and non-targeting or along the differentiation/pseudotime. It sounds like Fig3a shows only differentially expressed genes (between KO and control – are all of them used for the GSEA in 3b?) and how the expression of these genes changes (only in Gfra1 KO? cells) with pseudotime. Given the drastic changes in expression between pseudotime bins, which population was used for the differential expression testing between KO and control? All cells together? Only neurons?\nPage 10: “Collectively, these results support a role for the Gfra1 microexon in suppressing a transcriptomic program to prevent premature activation of neuritogenesis.” I am aware that also other microexon deletions have been shown to lead to a more neuron-like phenotype, however, this is still a bit counterintuitive to me. Forcing the non-neuronal splicing outcome (skipping) leads to a more neuron-like gene expression signature and phenotype. I would appreciate if the authors could elaborate on this in the discussion, especially given the fact that this seems to be a general consequence of neuronal nmicroexon deletion. In addition, an analysis of the commonalities and differences in the gene expression signatures of different microexon deletions would be helpful to strengthen the postulated convergence on FAK-Src pathway.\nPage 13: “A ‘module score’ representing expression of genes implicated in ASD…” – Are these genes implicated because of their premature expression (which would be in line with the transcriptomic signature of microexon deletion) or through loss-of-function/mutations? This distinction is relevant for the interpretation of the observed overlap in the list of genes.\nPage 13: “We observe that over pseudotime these premature maturation signatures show an early increase in expression of ASD-linked genes…” - Specifically ASD-linked genes? Or general neuronal and the observed association with ASD-linked genes is due to the fact that ASD-linked genes tend to have neuronal functions? What was the background gene list used in the enrichment analyses here?\nMinor comments:\nPage 3: “These challenges are particularly relevant to a program of evolutionarily conserved, 3-27 nucleotide-long ‘microexons’, which are predominantly spliced in the nervous system 14, 15.” – Should be “predominantly spliced in” or “predominantly included”.\nPage 13: “We also observe enrichment of ASD risk genes identified in three studies in the differential gene expression signatures detected in response to perturbation of Clasp1(A), Ptprf(A), Ptk2(A), Camta1 and Gfra1 microexons in our bulk RNA-seq” – should be clarified, e.g.: “…among the genes differentially expressed in response to…”\nPage 14: “Collectively, we observe that numerous neuronal microexons appear to function by ensuring robust and accurate timing of distinct stages of neurogenesis.” – Whether the word “numerous” is appropriate here is debatable. I would suggest to give specific numbers: How many microexons out of the 36 tested? How many microexons out of all known microexons?"}, "author_response": "Response to Reviewer Comments:\nReviewer #1\nIn this manuscript, Dupas and colleagues present CHyMErA-seq by extending the\nCHyMErA screening method developed earlier by the group. CHyMErA-seq allows one to\ndelete neural microexons in a high-throughput manner using CRISPR/paired gRNAs and\nevaluate the functional consequence of the genetic perturbation in individual cells by split-\npool single nuclei RNA-seq. This strategy is used to test a set of 37 conserved and neuron-\nspecific microexons in mESCs during in vitro neurogenesis. The study finds neuronal\ntranscriptomic shifts caused by microexon deletion and further examines the functional\nroles of two microexons (Gfra1 and Clasp1) in neurodevelopment. Deletion of microexons\nhas also been implicated in altering expression of ASD-risk gene networks, providing\ninsights into how dysregulated microexon splicing may contribute to disease. The\nfunctional evaluation of dynamic alternative splicing patterns, especially high-throughput\ntechnologies, is highly valuable. Indeed, within this class, the biological requirements of\nmicroexons are particularly important due to their potential connection with\nneurodevelopmental disorders.\nWhile the authors have made great strides in the execution of this high throughput strategy,\nthis reviewer has some concerns regarding the approach itself, and the robustness of the\ndownstream gene expression changes, and their relevance to maturation and disease. In\nparticular, the manuscript will benefit from strengthening key aspects, including microexon\ndeletions, as the entire body of work depends on this strategy, characterization of the\napproach, and phenotyping of neurons with exon deletions.\nWe thank the reviewer for their appreciation of the importance of the large-scale functional\nanalysis of microexons, interest in our strategy and results, and constructive comments,\nwhich we have addressed below.\nCHyMErA-seq screening:\n1. Microexon deletion: Figure 1, the authors show %del to account for deletions. This is\nlikely an overestimation as smaller PCR amplicon and typically favored. Even if one ignores\nthis bias, the estimated deletion frequency (e.g., top panels in Figure 1d) does not seems to\nbe consistent with the gel images in some cases.\nMoreover, the level of deletion in a bulk sample is not an especially useful metric, as what\nis more relevant is %cells with homozygous deletions. Carefully evaluating (for example\nwith ddPCR), and reporting these numbers will help the audience understand the\nrobustness of this dual Cas approach. The same problem applies to the application of\nCHyMErA-seq, where it is not clear to what levels homozygous microexon deletions\noccurs.\neXicient amplification of microexon-deleted amplicons, and that there are inherent\nchallenges in estimating editing eXiciencies in bulk cell populations, we would like to\nemphasize that the main goal of the experiments in Figure 1 was to optimize the CHyMErA\nsystem by assessing relative editing eXiciencies between diXerent Cas12a constructs.\nNevertheless, to address the reviewer’s first concern, percent editing eXiciencies indicated\nin Figure 1D do not appear congruent with apparent editing eXiciencies from visualizing the\ngel bands because a length correction was applied (see amendment to the Methods on\npage 21), which is summarized below) (Gonatopoulos-Pournatzis, Aregger et al. 2020). This\ncorrection accounts for weaker fluorescence signal due to reduced intercalation of\nethidium bromide into shorter amplicons. Notably, we also compared diXerences in\nestimated percent editing eXiciencies measured by PCR assays on bulk and ‘ground-truth’\nclonal cells using both corrected and uncorrected measurements (see Revision Fig. 1a\nbelow). We observe that the median absolute error after correction is lower (median =\ncorrected\n5.85, median = 9.22), thus validating the use of this approach.\nuncorrected\nPercent editing was calculated as the quotient of the length-normalized band intensities:\n𝐷𝑒𝑙𝑒𝑡𝑖𝑜𝑛𝐵𝑎𝑛𝑑\n!\"#$\"%!#&\n. 9\n𝐷𝑒𝑙𝑒𝑡𝑖𝑜𝑛𝐵𝑎𝑛𝑑\n'$\"(#)\n% Deletion = 100 *\n𝑈𝑛𝑐𝑢𝑡𝐵𝑎𝑛𝑑 𝐷𝑒𝑙𝑒𝑡𝑖𝑜𝑛𝐵𝑎𝑛𝑑\n!\"#$\"%!#& !\"#$\"%!#&\n. + 9\n𝑈𝑛𝑐𝑢𝑡𝐵𝑎𝑛𝑑 𝐷𝑒𝑙𝑒𝑡𝑖𝑜𝑛𝐵𝑎𝑛𝑑\n'$\"(#) '$\"(#)\nWe have also added text to the figure legend to highlight that this correction has been\napplied:\nFigure 1c and Extended Data Figure 1 a and 1b captions:\n“Editing e<iciency calculated from length-normalized quantification of ratios of edited and\ntotal band intensities (Methods)”\nTo address the reviewer’s second concern, we have assessed our method for estimating\nediting eXiciency from bulk genomic PCR by comparing measurements from ‘ground-truth’\nediting eXiciencies from genotyped single-cells derived from the bulk cell populations.\nImportantly, from these comparisons, we observe a significant correlation between bulk\nand ground truth editing eXiciencies (e.g. see Revision Fig. 1b below; R = 0.83, p = 0.012,\nPearson correlation), thereby supporting the validity of our approach.\ne=iciency. a) Comparison of the absolute di=erence in editing e=iciency measured using PCR\nassays performed on bulk compared to clonal cells, with (red) and without (blue) quantifications\ncorrecting for length of PCR products. b) Scatterplot comparing percent editing e=iciency\nestimated from genomic PCR assays performed on bulk cell populations, with editing e=iciency\nestimated from the proportions of edited cells observed from the analysis of the single cell clones\nderived from the CHyMErA deletions tested in this manuscript (red), and deletions used to generate\nadditional clonal lines (blue). R = 0.83, p = 0.012, Pearson correlation. N = 8 guides targeting 7\ndi=erent exons, median 12 clones per guide.\nWe thank the reviewer for recommending the use of ddPCR. While this technique could\npotentially provide more precise measurements of edited and unedited genomes, a\nstandard ddPCR assay would not be informative for distinguishing between homozygous\nand heterozygous deletions. This is because DNA is randomly partitioned in ddPCR\ndroplets and therefore lacks allelic linkage within single cells. As such, standard ddPCR\ncannot distinguish between, for example, a population of heterozygous cells and a 50:50\nmix of WT and homozygous-edited cells.\nTo address the reviewer’s comment, as mentioned above, we genotyped clonal cell lines\nderived from populations of cells edited with two independent guides with our enhanced\nCHyMErA system and observed that 37.5 and 50% of edited clones from these two\ndeletions harbour homozygous deletions. However, it should be noted that based on\nadditional analyses in our submitted manuscript and from the published literature (see\nbelow), we would expect that the majority if not all heterozygous deletions of microexons\nthat exhibit robust phenotypic changes when homozygously deleted to also result in\nmeasurable phenotypic changes. This conclusion is supported in the present study by the\nuse of our dCasRx-RBM25 targeting system for microexon modulation as an orthogonal\nassay to measure phenotypes. Specifically, partial skipping of microexons (eg. Clasp1(A),\nCamta1) directed by dCasRx-RBM25 yielded significant transcriptomic changes that\noverlap those detected in the CHyMErA-seq screen.\nresult in significant molecular and behavioural phenotypes. Particularly notable among\nthese are results from analyzing the eXects of partial loss of expression of the neuronal\nmicroexon-containing isoform of Cpeb4, which results in dysregulation of ASD genes and\nelectrophysiological and behavioural phenotypes in mice (Parras, Anta et al. 2018, Garcia-\nCabau, Bartomeu et al. 2024), and heterozygous germline deletion of a neuronal microexon\nin the Bak1 gene results in increased neuronal apoptosis during development (Lin, Zhang et\nal. 2020). Furthermore, in a previous study we demonstrated robust yet reduced\ntranscriptomic phenotypes when comparing heterozygous versus homozygous deletions of\na neural-diXerential exon in Puf60 in the same diXerentiation model as used in the present\nstudy (Han, Best et al. 2022).\n2. The representation of guide pairs in the sequenced nuclei library (and potential sources\nof bias) needs further illustration. The authors have targeted 37 exons with 9 guide pairs\neach. How many cells were recovered per guide pair for each target that passed QC should\nbe included. How is the correlation of DGEs between guide pairs targeting the same\nmicroexon aXected by cell#?\nWe refer the reviewer to Extended Data Figures 2a and 2b and the associated text (on page\n6) in which we describe the number of cells recovered for individual guides and target:\n“Each pair of microexon-targeting gRNAs was represented by a median of 44 cells, and\neach microexon target was represented by a median of 492 cells (Extended Data Figs. 2a,\nand 2b)”.\nTo address the impact of cell number on guide correlation, we plotted this relationship and\ncalculated the correlation (see Revision Fig. 2 below; R = 0.12, p = 0.35, Pearson\ncorrelation). This analysis reveals lack of a significant correlation between the numbers of\ncells expressing diXerent guides targeting the same microexon and changes in\ntranscriptomic signatures, thus indicating that cell number is not a driving factor in the\ncorrelations we are observing. To highlight this point, we have added the figure shown\nbelow to Extended Data Figure 4e and the following text (on page 7):\n“Importantly, cells with guides targeting the same gene exhibited significantly higher gene\nexpression correlation compared to both non-targeting (p = 7.50x10-12, Benjamini-Hochberg\ncorrected Wilcoxon rank sum test) and intergenic-targeting controls (p = 5.85x10-9,\nBenjamini-Hochberg corrected Wilcoxon rank sum test) (Extended Data Fig. 4d). This\ndegree of correlation is to a large extent independent of the number of cells recovered for\neach guide (Extended Data Fig. 4e).”\n4\n0.75\n0.50\nR = 0.12, p = 0.35\n0.25\n10 20 30 40 50 60\nNumber of cells per guide\nRevision Figure 2. Scatterplot comparing Pearson correlations for guides targeting the same\nmicroexon versus the number of cells represented by each guide. Pearson’s R and P value are\nshown.\n3. A major caveat of the CHyMErA-seq approach, when applied to mESCs, is the low editing\neXiciency, and cells with microexon deletion cannot be isolated directly. The use of QDA to\nclassify cells as KO vs non-perturbed based on gene expression profile is a creative way to\ndeal with the KO vs non-perturbed issue. However, this is an indirect workaround, and the\nreliability of the approach remains to be determined.\nWhile we acknowledge that a limitation of the CHyMErA-seq strategy is not being able to\ndirectly confirm microexon deletion in cells harbouring transcriptomic signatures, the\nvalidity of this screening approach is supported by the extensive follow-up experiments\nprovided in our study, including orthogonal assays performed on cell lines harbouring\nindependent, genotyped deletions (Figure 3), and results from independently perturbing\nmicroexons using dCasRx-RBM25 (Figure 4). We also note that estimating the ‘true’\nproportion of edited/eXected cells in bulk cell populations represents a more general\nchallenge encountered in CRISPR-based screens yet, as in the case of our work, has not\nhampered providing phenotypic read-outs that serve as a valuable basis of follow-up\nstudies (e.g. (Dixit, Pamas et al. 2016, Papalexi, Mimitou et al. 2021, Santinha, Klingler et al.\n2023, Jiang, Dalgarno et al. 2025) to list a few).\nThe authors compared CHyMErA-seq derived gene expression changes in response to\nmicroexon deletion and that determined by bulk RNA-seq of clonal samples, and the\ncorrelation is quite moderate (Figure 4d,e and Extended Data Figure 5d), which is very\nconcerning.\nWe understand the reviewer’s concern but given that we are comparing transcriptomic\nsignatures derived from single-nuclei RNA-seq and whole-cell bulk RNA-seq data, we\nwould not necessarily expect to observe strong gene-gene correlations. This is especially\nthe case when considering the considerable heterogeneity of the bulk cell populations and\n5\nnoitalerroc\nnosraeP\nalso the major technical diXerences involved, for example, in library-preparation and data\nprocessing methods. However, despite these significant diXerences, we observe positive\nand significant correlations between the bulk and single-nuclear RNA-seq profiles (Fig. 4d\nand 4e, Extended Data Fig. 5c). Additionally, we observe the same directional-shifts in the\nexpression of genes with a higher neuron diXerentiation rank in the bulk RNA-seq data (Fig.\n3c, Fig, 4b and 4c) as we do in the single nuclear signatures (Fig 2g). These findings\ntherefore support the overall consistency of the detection of altered gene expression\nprograms detected in both the snRNA-seq and the bulk RNA-seq data. This consistency is\nfurther reflected in a new analysis of enrichment of Gene Ontology terms associated with\nupregulated genes by directly comparing the single-nuclear and bulk RNA sequencing\ndata, which is now included as a new panel in Figure 3 (Figure 3b, shown also below).\nFor example, we observe enrichment of multiple common GO terms associated with neuron\ndiXerentiation and neuronal morphogenesis for genes up-regulated in both snRNA-seq and\nbulk RNA-seq after Gfra1 microexon deletion.\nPage 10:\n“Furthermore, we observe a significant overlap between the sets of Gene Ontology (GO)\nterms enriched for upregulated genes in response to Gfra1 microexon deletion in both bulk\nRNA-seq and CHyMErA-seq (odds ratio = 20.4, P = 2.2x10-16, Fisher’s exact test) (Fig. 3b).\nWe observe that most of the overlapping GO terms are related to processes involved in\ncellular morphogenesis and neuron di<erentiation.”\nhomophilic cell adhesion via plasma membrane adhesion molecules\nforebrain generation of neurons\ncell morphogenesis involved in neuron differentiation\neye development\nsensory organ development Method\nregulation of nervous system development\nsystem development bulk\nmulticellular organism development\nlocomotion snRNA−seq\ncell development\nregulation of multicellular organismal process\nplasma membrane bounded cell projection organization\np < 0.05\ncell migration\nmulticellular organismal process\np < 0.01\nanatomical structure development\ndevelopmental process\ncell communication\nsignaling\nsignal transduction\nlog(odds ratio)\n1 2 3\nRevision Figure 3. Heatmap of log-transformed odds ratios of Gene Ontology (GO) biological\nprocess terms significantly enriched (P < 0.05, Fisher’s exact test) among up-regulated genes\nidentified by bulk RNA-seq and CHyMErA-seq following Gfra1 microexon deletion. P values were\ncorrected by Benjamini-Hochberg correction.\n6\n4. It seems that no replicate screening was performed, which adds additional challenge to\nevaluate the robustness of the screening.\nWe performed the CHyMErA-seq screen using multiple distinct hgRNAs targeting each\nmicroexon to serve as replicates for edits and to enable assessment of the robustness our\ntranscriptomic measurements. This approach for assessing replication of measurements\nhas the major advantage of avoiding potentially confounding variability due to possible\nbatch eXects that one would expect between replicates of entire screen experiments that\ninvolve a large number of procedures and months of work-up. We also emphasize again\nthat our study incorporates extensive validation through the application of dCasRx-RBM25\nmicroexon perturbations and other follow-up assays, which served to both confirm and\nextend results from the CHyMErA-seq data.\n5. In some parts, the description of bioinformatics analysis is not very clear. For example, it\nis diXicult to figure out how perturbation calling was done precisely. The authors might\nconsider using less narratives and include some mathematical descriptions.\nWe thank the reviewer for their comment and have updated the Methods to provide a more\ndetailed explanation of perturbation calling. In addition to including equations describing\nediting eXiciency calculations, we have provided the following new text and associated\nmathematical equations that describe the perturbation calling strategy:\nPage 27:\n“We calculated a normalized gene expression vector for each cell (𝑥@ ) subtracting the\n!\naverage expression of the neighboring cells B𝑥̅ D for each gene:\n*\n!\n𝑥@ = 𝑥 −𝑥̅\n! ! *\n!\nwhere:\n∑ 𝑥\n𝑥̅ =\n+∈*! +\n* ! 𝑘\nPage 28:\n“To classify cells as perturbed or non-targeting we applied quadratic discriminant analysis\n(QDA), implemented in sci-kit learn (v1.2.2), to calculate classification assignments from\nmultivariate Gaussians fit on normalized expression values 𝑥@ . QDA models each class as a\n!\nmultivariate Gaussian with its own class-specific covariance matrix Σ and mean µ and\n- -\nassigns classifications using Bayes’ rule:\n7\n𝑃(𝑥@ ∣∣ 𝑦 = 𝑘)⋅𝑃(𝑦 = 𝑘)\n! ! !\nP(𝑦 = 𝑘 ∣∣ 𝑥@ ) =\n! ! 𝑃(𝑥@ )\n!\nfor which we calculate the class conditional likelihood as:\n1 1\nP(𝑥@ ∣∣ 𝑦 = 𝑘) = expU− (𝑥@ −µ )2Σ31(𝑥@ −µ )V\n! ! (2π)./0|Σ |1/0 2 - - -\n-\n“For each cell, an empirical p-value was calculated by comparing its observed classification\naccuracy to the null distribution.\n1+∑* 𝐼Y𝑝̂ ≤ 𝑝null,\"\n𝑝 = \";1 ! !\n45'6$ ! 𝑁+1\nwhere 𝑝̂ is the observed classification accuracy of cell 𝑖, and 𝑝 \"6'',\" is the accuracy in the\n! !\nnth permutation.”\nPhenotyping:\n6. Are the 37 microexons frame preserving? This reviewer assumes so, but it will be helpful\nto clearly indicate this.\nWe have updated the text on page 5 to clarify:\n“All of the targeted microexons are frame-preserving and expected to insert 3-9 amino\nacids into the coding sequence”\n7. The authors have positive control gRNA pairs that target constitutive exons, but do not\nprovide much detail on the transcriptomic characterization of these cells.\nInclusion of the constitutive exon-targeting guides was intended to serve as positive\ncontrols to benchmark transcriptomic responses relative to microexon deletion, intergenic\nand intronic guide controls. However, assessing transcriptomics signatures from these\nguides alone was challenging due to the relatively small numbers of cells recovered from\nthese deletions, which was in part due to constraints imposed by a size limitation of our\nhgRNA library design, as we prioritized the library for targeting microexons over constitutive\nexons in the same genes. However, in response to the reviewer’s comment, we have\n8\nperformed additional analyses to assess the extent to which transcriptomic signatures\ndetected following microexon deletion are detected in cells expressing guides targeting\nconstitutive exons (‘exonic’ guides) designed to generate gene knockouts.\nFor each microexon deletion that resulted in a detectable transcriptomic phenotype, we\nidentified diXerentially expressed genes (FDR < 0.1) and performed PCA on their\nnormalized expression across single cells (see ‘Normalization to non-targeting cells’ in the\nMethods). We then compared transcriptomic signatures of microexon-targeting, exonic or\nintergenic negative control guides, with those from non-targeting cells. We calculated\nHotelling’s T² statistic on the top 40 principal components, or fewer if limited by the\nnumber of cells in each group, to assess statistical significance.\nThis analysis reveals that, among all of the targeted genes, no constitutive exon-targeting\nguides elicited significantly altered the same transcriptomic signatures that were detected\nin cells targeted for microexon-deletion (i.e. relative to non-targeted cells) (Revision Figure\n4, see below). A similar result was observed when comparing with intergenic and non-\ntargeting negative control cells. These results suggest that the transcriptomic signatures\ninduced by microexon deletion are distinct from those generated by constitutive exon-\ntargeting guides, consistent with splicing dysregulation-specific phenotypes. We have\nadded the following text to the manuscript and the accompanying Extended Data Fig. 4g to\naddress the reviewer’s comment.\nPage 7:\n“To assess the specificity of microexon-deletion phenotypes, we performed pairwise\ncomparisons of microexon deletion-associated transcriptomic signatures between cells\nexpressing microexon-targeting hgRNAs or cells expressing guides targeting constitutive\nexons in the same gene, with cells expressing non-targeting guides (Methods) (Extended\nData Fig. 4g). This analysis reveals that transcriptomic signatures derived from microexon-\ndeletions are distinct from those detected in response to guide-directed perturbation of the\ncorresponding genes.”\nexpressing control guides or guides targeting constitutive exons in the same genes. Heatmap\nshowing -log10(corrected p-values) for pairwise comparisons of transcriptomes between non-\ntargeting cells and microexon-deletion cells (targeting), constitutive exon-targeting (exonic) or\nintergenic negative control guides (grey indicates NA). Cells were embedded using principal\ncomponent analysis of di=erentially expressed genes and Hotelling’s T2 test was used to assess\ngroup-level di=erences in multivariate expression space. P* > 0.05, P** > 0.01, P*** > 0.001;\nBenjamini-Hochberg corrected Hotelling’s T2 test.\n8. Pseudotime analyses indicates that the KO cells are broadly enriched in the\nneuroepithelia clusters (Fig 2d). This suggests there are diXerentiation defects and KO cells\nare stuck in the neuroepithelium state? As the microexons are typically not included to\nsignificant levels at this stage of diXerentiation, how do the authors explain KO cells\nenrichment in neuroepithelia. Moreover, this seems incongruent with the idea that\nmicroexons deletions increase maturation.\nIn fact, the analyzed microexons are strongly included at this stage (e.g. refer to Day In Vitro\n7 timepoint in Figure 1f, shown below). The pseudotime analysis merely indicates that\nwithin the trajectory these cells are less mature. Contrary to the reviewer’s interpretation\nthat these cells may be arrested at this stage of diXerentiation, these cells have an\nupregulated diXerentiation signature at this early stage and do diXerentiate. This is\nillustrated in the validation experiments focusing on the Gfra1 microexon deletion,\nin which we observe that the microexon deletion does not cause diXerentiation arrest, but\nrather an upregulation of diXerentiation-related transcriptional programs and increased\nneuritogenesis (refer to Figure 3).\nseq analysis during in vitro neuronal di=erentiation. From Figure 1f.\n9. More can be done to solidify/validate gene expression changes observed in the snRNA\ndata. While the authors decided to focus on Gfra1 and Clasp1 for additional analyses and\nmade eXorts to generate clonal cell lines with microexon deletion (or skipping/inclusion\nusing dCasRx-RBM25). The evidence provided to support the functional role of each exon is\nless than compelling in the current form. For example, for Gfra1, the transcriptomic\nchanges inferred from the screen essentially has no correlation with that from analysis of\nclonal cell lines (Extended Data Figure 5d). For Clasp1, the impact of microexon on\ndownstream phosphorylation appears to be quite moderate.\nTo validate CHyMErA-seq-detected gene expression changes, we performed bulk RNA-seq\nanalysis on clonal cell lines harbouring microexon deletions that were generated using\ngRNAs that are diXerent from those used for the CHyMErA-seq screen, and we also\nmanipulated microexons using dCasRx-RBM25 targeting. As described in response to the\ncomment # 3 above, gene-gene level correlations are not expected to be high between\nsnRNA-seq and bulk RNA-seq data, in large part because of the significant technical\ndiXerences involved in sample preparation and analysis. Nevertheless, we still observe\n11\nsignificant and positive correlations (see Extended Data Figure 5c, Figure 4d and 4e).\nFurthermore, consistently, we also observe significant overlaps between GO functional\nterms that are enriched in the sets of genes with neuronal diXerentiation-related\ntranscriptional signatures between these datasets (refer to Revision Figure 3). By\nperforming the same analysis on gene expression signatures resulting from Clasp1\ndCasRx-RBM25 manipulation, we also observed significant overlap with GO term\nenrichments derived from CHyMErA-seq and the bulk RNA-seq (odds ratio = 52.1,\nNeurons\nP < 2.2x10-16; odds ratio = 29.3, P < 2.2x10-16, Fisher’s exact test) (see Revision\nNeurons mESC mESCs\nFig.6 below).\nPage 11:\n“We further observe significant overlap of GO terms enriched among genes that are\nupregulated in the CHyMErA-seq gene signature, and in mESCs (odds ratio = 29.3, P <\n2.2x10-16, Fisher’s exact test) and DIV7 neurons (odds ratio = 52.1, P < 2.2x10-16, Fisher’s\nexact test)) (Fig. 4f). These GO terms are associated with functions related to the regulation\nof nervous system development, morphogenesis and cell-signaling.”\nTimepoint\nbrain development\ncentral nervous system development\nNeurons\nanimal organ morphogenesis\ncell adhesion\nmESC\nsystem development\nregulation of developmental process\nsnRNA−seq\ncell development\ncellular developmental process\npositive regulation of multicellular organismal process\nanatomical structure morphogenesis Method\nmulticellular organism development\nanimal organ development bulk\nmulticellular organismal process\ncell−cell signaling snRNA−seq\npositive regulation of transcription by RNA polymerase II\ncell communication\nsignaling\nregulation of signaling p < 0.05\nregulation of cell communication\nsignal transduction p < 0.01\nlog(odds ratio)\n0.5 1 1.5\nRevision Figure 6. Heatmap of the log(odds ratio) of representative significantly enriched GO terms\nin bulk RNA-seq and CHyMErA-seq following Clasp1(A) microexon exclusion compared to\ninclusion in mESCs and neurons (Methods). P values are FDR adjusted.\n10. The authors used a KS-test to evaluate alterations of gene expression correlation with\nneuron diXerentiation rank. It is quite diXicult to understand what these significant\ndiXerences really mean.\nof a gene correlates with the position of a given cell on the diXerentiation trajectory. These\nhighly ranked genes are enriched in neuronal processes, while lowly ranked genes are\nenriched in processes involved in cell division and catabolism (see Extended Data Fig. 4i).\nWe agree with the reviewer that the use of the KS-test here is somewhat indirect and have\nadapted the analysis with a more straightforward comparison. We now assess the\nsignificance in the distribution of neuron diXerentiation ranks for up- and down-regulated\ngenes with a two-sided Wilcoxon rank sum test comparing the up- or down-regulated\ndistribution of ranks with the distribution derived from genes that are not significantly\naltered (i.e. FDR ³ 0.05). Therefore, we can interpret the enrichment of genes based on their\nregulation towards higher or lower neuron diXerentiation rank scores. We have added the\nfollowing statement to the text to clarify this analysis and have accordingly updated Figures\n2g, 3c, 4b and 4c to reflect the new statistics:\nPage 9:\n“To assess how microexon deletion eXects neurogenesis, we compared the distributions of\nscores for up- and down-regulated genes for each perturbation, or in the case of Ptk2,\npooled scores for genes a<ected across the three microexon perturbations, due to the\nrelatively low cell numbers and sparse data recovered for each individual Ptk2\nperturbation. Notably, compared to background genes (i.e. non-di<erentially expressed,\nFDR ³ 0.05), we observe that upregulated genes associated with altered transcriptomic\nsignatures are significantly enriched for higher CellRank2 scores for microexon deletions in\nthe Bin1, Clasp1, Gfra1, Med23, Ptprf and Ralgapb genes (P < 0.05, Benjamini-Hochberg\ncorrected two-sided Wilcoxon rank sum test) (Fig. 2g).”\n11. Ptk2 microexons are known to dramatically elevate autophosphorylation. Do the\nauthors detect signatures of modulated Ras signalling/other downstream signaling activity?\nSimilarly, do the authors observe elevated morphology or migration changes downstream of\nPtk2 microexon modulation?\nWe are aware of the excellent work performed on the Ptk2 microexons and have referenced\nthese papers (Burgaya, Toutant et al. 1997, Toutant, Costa et al. 2002). We refer the\nreviewer to page 13 where we highlight the link between microexon deletions, including\nPtk2 microexons, and the Ras/MAPK-signaling pathways (which is downstream of FAK-\nintegrin signaling):\nPage 13:\n“Of these six dysregulated regulons associated with microexon deletion, three are\ndownstream of MAPK-ERK signaling pathway (Ets2, Erg1 and Elk1), which has been\nimplicated in core processes during neurogenesis (Gantner, Hunt et al. 2021). Notably, in\naddition to the role of the Clasp1(A) microexon in augmenting MAPK-ERK signaling shown in\nthis study, roles in MAPK signaling have been previously described for Bin1 (Wang, Jia et al.\n2017), Gfra1 (Paratcha, Ledda et al. 2003), Med23 (Wang, Huang et al. 2009, Yang, Zhao et\n13\nal. 2012, Yang, Li et al. 2023), Ptprf (Wang, Weng et al. 2000, Wang, Kuo et al. 2007) and\nPtk2 (Mitra, Hanson et al. 2005) proteins (which also form a protein-protein interaction\nnetwork, Fig. 2f) suggesting that the microexons in these proteins may convergently\nregulate this signalling pathway.”\nWe observe significant downregulation of the regulon score associated with Egr1, which is\nactivated downstream of MAPK signaling pathways in response to Ptk2 microexon\ndeletions (Figure 5a). Since Ptk2 microexons have already been characterized in some\ndetail, we did not follow up on these with morphological assays and instead focused on\nrelatively less-well characterized microexons in the Gfra1 and Clasp1 genes.\n12. GDNF signaling is known to promote survival and growth of neurons. The authors find\nthat Gfra1 microexon deletion modestly increases neurite growth and stimulates\nexpression of mature neuronal genes. Is this a consequence of increased GDNF signalling?\ni.e., does microexon inclusion decrease Gfra1 signaling? Can the authors supplement\nneurons with GDNF to amplify the neuritogenesis phenotypes?\nWe interpret that these changes are due to increased Gfra1 signaling due to previous work\nshowing that inclusion of the microexon in Gfra1 decreases its aXinity for Gdnf and Nrtn\nand thereby downregulates Ret-signaling (Charlet-Berguerand, Le Hir et al. 2004). We\ninadvertently omitted this citation from the original text and have corrected this error in the\nrevised manuscript. While we agree with the reviewer that it would be interesting to\ntitrate\nGdnf and measure eXects on neurite outgrowth, we feel that this experiment is beyond the\nscope of our manuscript, which is primarily focused on introducing the CHyMErA-seq\ntechnology, validation of this system, and initial follow-up characterization on selected\ninteresting screen hits.\n13. Camta1 microexon does not appear to be included in mESCs. Nevertheless, the\nchanges in gene expression upon deletion in mESCs appears to enrich downregulated\ngenes that are aXected in autism (Fig 5d). this feature appears significant only in the\nmESCs and not in neurons. Similar results are also obtained for Clasp1 microexon as well.\nHow do the authors explain observations of transcriptomic phenotypes in mESCs that do\nnot express much of the microexon, but not in neurons from the same founders?\nActually, the Camta1 and Clasp1 microexons do show basal inclusion in mESCs\n(approximately 56 PSI and 31 PSI, respectively; see Revision Figure 7 below, and Revision\nFigure 5 (i.e. Fig 1f)). Furthermore, in these experiments we are comparing dCasRx-RBM25-\ndirected skipping and inclusion in these cells, which aXords a larger diXerence in the\ninclusion levels of these microexons between the two conditions. Strikingly, we observe\nthat dCasRx-RBM25-directed skipping of the microexon in mESCs is suXicient to induce\ntranscriptomic changes corresponding to a more neuronal signature, as compared to\ndCasRx-RBM25-directed exon-inclusion mESCs. These results highlight that by forcing\nsplicing level diXerences of these microexons in mESCs using dCasRx-RBM25, we observe\ntranscriptomic changes similar to those seen in the CHyMErA-seq screen data, which are\n14\nenriched in ASD-linked genes that are often implicated in early nervous system\ndevelopment (Satterstrom, Kosmicki et al. 2020) as is indicated in Figure 5d.\nRevision Figure 7. RT-PCR quantification of inclusion levels of Camta1 and Clasp1(A) microexons\nwhen transduced with non-targeting (‘NT’) exon skipping (‘EX’) or inclusion (‘INC’) guides in\ndCasRx-RBM25 expressing mESCs.\nOther minor points:\n14. Explain hgRNA\nWe have elaborated on the description of hgRNAs in the Introduction:\nPage 3:\n“CHyMErA employs simultaneous expression of Cas9 and Cas12a nucleases, together\nwith hybrid guide (hg)RNAs, which are fusions of Cas9 and Cas12a guides expressed from a\nsingle U6 promoter and processed into pairs of individual Cas9 and Cas12a gRNAs that\ndirect genomic fragment deletion.”\n15. Line 689, missing “)” after CRISPOR.\nCorrected.\n16. Line 755, paring --> pairing\nCorrected.\n17. Regulon scores were z-normalized and p values **and** were calculated for each\npairwise comparison by two-tailed unpaired t-test followed by Benjamini-Hochberg\ncorrection.\nCorrected.\n18. Fig. 2b, colors for some of the clusters are too light to be visible\n19. Fig. 2e, what are the columns?\nWe have replaced this figure with a revised version that illustrates enrichment of gene\nontology terms for the DEGs associated with microexon deletions thereby providing a more\nbiologically interpretable figure.\n20. Fig. 4 legend, “Empirical cumulative distribution function (ECDF) of diXerentially\nexpressed genes (DEGs, FDR < 0.05) associated with microexon inclusion compared to\n**inclusion** for mESCs.”\nCorrection made.\nReviewer #2\nThis manuscript presents CHyMErA-seq, a platform designed to profile transcriptomic\nchanges at single-cell resolution in response to targeted splice isoform perturbations.\nApplying this approach to 36 neuronal microexons, the authors identify consequences of\nthese alternative splicing events on gene expression dynamics during neurogenesis.\nSeveral microexon-containing genes tested here are implicated in neuron projection\ndevelopment and converge on signaling pathways such as FAK-Src. Functional analyses\nreveal that specific microexon deletions studied here (similar to the eXect of microexons\nstudied elsewhere) lead to increased neurite outgrowth and in general a premature\nactivation of neuronal gene expression programs. The authors also observe parallels\nbetween microexon-dependent transcriptomic signatures and gene expression patterns\nobserved in autism spectrum disorder, in line with previous studies.\nThis study builds on a large body of research from the Blencowe and other labs, which\nprovided the basis for the experimental tools applied here. In the present study, the authors\nimproved editing eXiciency and coupled the CHyMErA pipeline with scRNA-seq readout. In\naddition, most if not all the main results obtained using scRNA-seq here confirm or are\nin line with previous studies, as the authors themselves point out repeatedly (activation of\nneuronal gene expression signatures, increased neurite outgrowth, link to ASD,…). On the\none hand, this is an encouraging proof-of-concept of the experimental approach, but on\nthe other hand raises questions about the novelty of the findings. In my opinion, it will be\nnecessary to delineate more specifically which results are novel, which are confirmatory,\nand where the single-cell and isoform-level perspective (which is the main “selling point” of\nthis study) allowed new biological insights into a previously studied biological question.\nSpecific points and suggested analyses are detailed below. Especially considering the\nnecessity of pointing out the novelty of the findings compared to previous studies, it would\n16\nbe advisable to be more specific in the wording of the results and avoid vague statements\nsuch as the ones below:\nPage 12: “Overall, these data confirm transcriptomic signatures detected in the CHyMErA-\nseq screen and further support important new functions of microexons in neurogenesis”\nPage 14: “Collectively, these findings highlight the importance of isoform-specific protein\nfunctions of genes in the control of critical and dynamic processes during neurogenesis,\nand their dysregulation in neurodevelopmental disorders.”\nPage 15: “Applying this system to neuronal microexons revealed unexpected functional\nproperties [which ones?] of this important class of alternative splicing in the context of\nneurogenesis, as well as new insights into the multilayered links between microexons and\nbrain disorders [I do not agree that the manuscript contains enough novel insights in this\nrespect to justify this sentence].”\nWe thank the reviewer for their feedback and recommendation to better highlight the novel\nand distinct aspects of our study. In our view, one of the most important and novel\ncontributions of our study is that it reports the first scalable system for the systematic\nperturbation of exons coupled to a single cell transcriptomics read-out. This approach\nenables for the first time the parallel interrogation of the functions of dozens of exons in a\nsingle experiment, which in turn aXords direct comparisons of the activities of multiple\nsplice variants in specific biological processes and pathways. In particular, while a small\nnumber of previous studies have revealed functions of microexons in the control of\ntranscriptomic outputs and neuritogenesis, these studies focused on analyzing one\nmicroexon at a time in bulk cell populations and at single time points. In contrast, by\nanalyzing the functions of multiple microexons at single cell resolution at diXerent stages\nof development our study reveals convergent functions of microexons in specific signalling\npathways and in the timing of neurogenesis-associated transcriptomic signatures. Lastly,\nwhile our study purposefully targeted a few microexons that have been previously studied\nto compare specific biological outcomes (i.e. aXording a level of independent validation),\nmost of the microexons interrogated in our study have not been previously analyzed on any\nlevel and therefore their functional roles were entirely unknown.\nTo address the reviewer’s comments relating to sentences on pages 12, 14 and 15, we have\nmodified the text to clarify the most important and novel findings or out study, as well as\nhow our results relate to – but complement and extend – the findings of previous studies:\nPage 12:\n“Overall, these data confirm the dysregulation of gene programs associated with neuron\nmorphology and development observed in the CHyMErA-seq screen and extend these\nobservations by highlighting disruption of downstream signaling cascades as a result of\nmicroexon perturbation.”\n“Taken together with the observations described above indicating functions for these\nmicroexons in controlling neuronal development, our results highlight a potential role for\ntheir dysregulation in contributing to the etiology of neurodevelopmental disorders.”\nPage 16:\n“Applying this system to neuronal microexons uncovered alternative splicing-dependent\ngene expression signatures associated with the regulation of developmental timing during\nneurogenesis, and further illuminated an additional connection between microexon\ndysregulation and neurodevelopmental disorders.”\nFurthermore, to provide context and contrast our work to previous studies, we have added\nthe following text to the discussion:\nPage 16:\n“In contrast to these previous studies, through the development of the CHyMErA-seq\nsystem the results of the present study revealed new and unexpected functions of\nmicroexons in the control of single cell transcriptomes during neurogenesis.”\nSpecific comments:\nPage 8/Fig 2: “The Ptk2 microexons (denoted as ‘A’, ‘B’ and ‘C’) display coordinated splicing\npatterns in vivo…”: Their splicing is coordinated, but here the authors delete one of them at\na time: Consequently, deletion of microexon alone will create an artificial protein that is not\nusually produced (because they are included or excluded together). Is there in vivo\nevidence for the existence of transcripts with two of the three microexons included and one\nskipped? In addition, the authors state that “notably, we observe transcriptomic signatures\nfor all three microexons targeted in the Ptk2 gene (Fig. 2e)”. First of all, this should be\nrephrased to be less vague and more accurate: (a) …transcriptomic signatures of the\nDELETION of all three microexons…; (b) …we observe DISTINCT transcriptomic\nsignatures…; (c) these apparent diXerences between the eXect of the deletion of the three\nmicroexons should be discussed (in particular, B and C are located in the same domain,\nbut yet deletion of each of them shows very diXerent eXects on the transcriptome).\nWe thank the reviewer for these important questions. While there are predicted transcripts\nin Refseq for Ptk2 isoforms matching all possible combinations of alternative splicing\nexons A, B and C, some of these transcripts appear to represent minor fractions of the total\ntranscript pool for this gene. This is especially true of transcripts containing exons B or C,\nwhich show coordinated splicing (‘A’-only transcripts are annotated in GENCODE and show\nan earlier inclusion in neurodevelopment). To directly observe the presence of B and C\ncombinations in neural data, we analysed single-cell isoform sequencing data from a\n18\nrecent paper that utilized samples from the developing human brain (Patowary, Zhang et al.\n2024) (Revision Figure 8 below). This analysis indicates that 20% of transcripts containing\n‘B’ or ‘C’ harbour one or the other microexon, while 80% harbour both. However, given the\nframe-preserving lengths of these exons (18 and 21 nucleotides), we expect that the\nresulting transcripts are still translated into proteins.\nRevision Figure 8. Bar plot of Ptk2 exons present in long-read single-cell data from developing\nhuman embryos (Patowary, Zhang et al. 2024). Bars represent the percent of transcripts containing\nPtk2 microexons ‘B’ and/or ‘C’ (n = 60 total transcripts).\nFinally, while we individually targeted the Ptk2 A-C microexons, there were relatively few\ncells harbouring these deletions resulting in sparse data for the individual perturbations.\nThe description of the transcriptomics read-out therefore relates to the average changes in\ngene expression after pooling data from cells harboring the diXerent perturbations. We\nhave modified the text to clarify this:\nPage 9:\n“To assess how microexon deletion eXects neurogenesis, we compared the distributions of\nscores for up- and down-regulated genes for each perturbation or in the case of Ptk2,\npooled scores for genes a<ected across the three microexon perturbations, due to the\nrelatively low cell numbers and sparse data recovered for each individual Ptk2\nperturbation.”\nme if the focus here is on gene expression changes between Gfra1 microexon knockout and\nnon-targeting or along the diXerentiation/pseudotime. It sounds like Fig3a shows only\n20\ndiXerentially expressed genes (between KO and control – are all of them used for the GSEA\nin 3b?) and how the expression of these genes changes (only in Gfra1 KO? cells) with\npseudotime. Given the drastic changes in expression between pseudotime bins, which\npopulation was used for the diXerential expression testing between KO and control? All\ncells together? Only neurons?\nTo address the reviewer’s question, we have replaced the previous Figure 3a with a simpler\nversion that is also shown below (Revision Figure 9). The figure now shows the scaled\naggregate expression of DEGs derived from Gfra1 microexon deletion in the CHyMErA-seq\nscreen over pseudotime for both Gfra1 deletion cells and non-targeting cells. We hope that\nthis panel better illustrates the temporal changes in the expression of upregulated genes as\na consequence of the Gfra1 microexon deletion cells through pseudotime, as compared to\nnon-targeting cells. The GSEA was performed using log-fold changes between Gfra1\ntargeted cells and non-targeting cells. We have moved this figure to Extended Data Figure\n5b. The legends for these figures (and others in our manuscript) have been expanded in\norder to make them easier to follow and the corresponding procedures have been added to\nthe Methods:\n“Temporal analysis of di1erential gene expression\nCell-wise scores computed per cell by summing normalized log-expression values (log1p-\ntransformed) of genes significantly upregulated in Gfra1 microexon deletion (KO) versus\nnon-targeting (NT) cells (logFC > 0, FDR < 0.05). Scores were z-scaled across all cells.\nA generalized additive model (GAM) was then fitted to the scaled scores as a function of\npseudotime using a smooth interaction with genotype. Model fits and 95% confidence\nintervals were computed using the mgcv package in R(Wood 2011).”\nRevision Figure 9. Temporal dynamics of Gfra1 microexon deletion-associated transcriptomic\nsignature. A generalized additive model (GAM) was fitted to scaled expression scores for genes\nsignificantly upregulated (logFC > 0, FDR < 0.05) in Gfra1 microexon deletion cells (KO, n = 103; red)\nand non-targeting controls (NT, n = 278; grey) across pseudotime (adjusted R² = 0.72; KO vs. NT\nestimate = 0.39; p = 6.7 × 10⁻⁹) (Methods). Shaded areas represent the 95% confidence intervals of\nthe fit.\ntranscriptomic program to prevent premature activation of neuritogenesis.” I am aware that\nalso other microexon deletions have been shown to lead to a more neuron-like phenotype,\nhowever, this is still a bit counterintuitive to me. Forcing the non-neuronal splicing outcome\n(skipping) leads to a more neuron-like gene expression signature and phenotype. I would\nappreciate if the authors could elaborate on this in the discussion, especially given the fact\nthat this seems to be a general consequence of neuronal nmicroexon deletion.\nWe have elaborated on these findings in the Discussion by providing the following text:\nPage 15:\n“Collectively, we observe that numerous neuronal microexons appear to function by\nensuring robust and accurate timing of distinct stages of neurogenesis. In particular,\nseveral of the analyzed microexons have convergent e<ects on transcriptomes that suggest\nthey serve as molecular ‘brakes’ during neurogenesis. This unexpected observation\nsupports the conclusion that specific microexons may have important roles in modulating\nprotein functions to meet context-specific requirements during neurogenesis. For example,\nthe Gfra1 microexon reduces ligand-binding a<inity(Charlet-Berguerand, Le Hir et al. 2004),\nand its deletion in this study leads to the increased expression of genes involved in\nneuronal di<erentiation and morphogenesis, along with enhanced neuritogenesis. In\ncontrast, deletion of Ptk2 microexons (which normally enhance phosphorylation\nactivity(Burgaya, Toutant et al. 1997, Toutant, Costa et al. 2002)) leads to downregulation of\nneurogenesis-related genes (Fig. 2g). This is opposite to the e<ect seen with deletion of the\nClasp1 microexon, which promotes FAK phosphorylation and accelerates neurogenesis.”\nIn addition, an analysis of the commonalities and diXerences in the gene expression\nsignatures of diXerent microexon deletions would be helpful to strengthen the postulated\nconvergence on FAK-Src pathway.\nWe have updated Figure 2e (Revision Figure 10, below) to highlight the functional terms\nenriched for each microexon deletion gene expression signature in an eXort to enhance the\ninterpretability of these transcriptomic signatures. This analysis shows that several of the\nmicroexon deletions result in dysregulation of genes involved in gene sets broadly involved\nin nervous system development and neuron projection guidance.\nRevision Figure 10. Heatmap of top 10 (sorted by FDR adjusted p-value) gene ontology (GO,\nbiological processes) terms for di=erentially expressed genes (FDR < 0.05) resulting from microexon\ndeletions. GO terms are grouped according to similarity and colored according to parent GO terms.\nOutline indicates Benjamini-Hochberg corrected P < 0.05. Grey indicates ‘NA’.\nPage 13: “A ‘module score’ representing expression of genes implicated in ASD…” – Are\nthese genes implicated because of their premature expression (which would be in line with\nthe transcriptomic signature of microexon deletion) or through loss-of-function/mutations?\nThis distinction is relevant for the interpretation of the observed overlap in the list of genes.\nThese genes are not implicated in ASD based on expression, but rather their genetic link to\nASD. The point we are making here is that ASD-linked genes fall within this program of\ndysregulated genes upon microexon deletion. This is clarified on page 13:\n“A ‘module score’ representing expression of genes with genetic links to ASD (using ASD\nassociation evidence scores 1-3 from Simons Foundation Autism Research\nInitiative (SFARI) (Abrahams, Arking et al. 2013))”.\nPage 13: “We observe that over pseudotime these premature maturation signatures show\nan early increase in expression of ASD-linked genes…” - Specifically ASD-linked genes? Or\ngeneral neuronal and the observed association with ASD-linked genes is due to the fact\nthat ASD-linked genes tend to have neuronal functions? What was the background gene list\nused in the enrichment analyses here?\nASD-linked genes represent a subset of the dysregulated genes with neuronal functions.\nWe have added the following text to clarify this point:\n22\nnoiteled\nnoxeorciM\ncell morphogenesis involved in neuron differentiation epithelial tube morphogenesis\npositive regulation of synaptic transmission cell growth\nregulation of nervous system development synapse assembly\ncell−cell adhesion developmental pigmentation\nsensory organ development other\ncell recognition\nGfra1\nRalgapb\nPtprf\nBin1\nClasp1\nMed23\nCamta1\nPtk2\n−log10(FDR)\n0 2 4 6\nPage 14:\n“We observe that over pseudotime these premature maturation signatures include an early\nincrease in expression of ASD-linked genes (pseudotime bins 3-6. p < 0.05, Benjamini-\nHochberg corrected two-sided Wilcoxon rank sum test) followed by a lagging increase in\nnon-targeting cells (pseudotime bin 7, p < 0.05, Benjamini-Hochberg corrected Wilcoxon\nrank sum test), before reaching a plateau (Fig. 5c).”\nThe background set in the enrichment analysis was all expressed genes used in the\ndiXerential expression analysis.\nMinor comments:\nPage 3: “These challenges are particularly relevant to a program of evolutionarily\nconserved, 3-27 nucleotide-long ‘microexons’, which are predominantly spliced in the\nnervous system 14, 15.” – Should be “predominantly spliced in” or\n“predominantly included”.\nWe have corrected this to:\nPage 3:\n“These challenges are particularly relevant to a program of evolutionarily conserved, 3-27\nnucleotide-long ‘microexons’, which are predominantly included in the nervous system”.\nPage 13: “We also observe enrichment of ASD risk genes identified in three studies in the\ndiXerential gene expression signatures detected in response to perturbation of Clasp1(A),\nPtprf(A), Ptk2(A), Camta1 and Gfra1 microexons in our bulk RNA-seq” – should be clarified,\ne.g.: “…among the genes diXerentially expressed in response to…”\nWe have made this correction\nPage 14: “Collectively, we observe that numerous neuronal microexons appear to function\nby ensuring robust and accurate timing of distinct stages of neurogenesis.” – Whether the\nword “numerous” is appropriate here is debatable. I would suggest to give specific\nnumbers: How many microexons out of the 36 tested? How many microexons out of all\nknown microexons?\nWe have changed ‘numerous’ to ‘multiple’. However, since 10 out of 36 of the microexon\ndeletions tested that showed robust transcriptomic changes in our analysis displayed\naltered timing of neurogenesis-associated gene expression patterns, we consider it likely\nthat this functional property of microexons may be more broadly relevant. We have\nmodified the text on page 15 to reflect this distinct observation relating to the function of\nmicroexons in our study.\nPage 15:\n23\n“Collectively, we observe that multiple neuronal microexons appear to function by ensuring\nrobust and accurate timing of distinct stages of neurogenesis…”\nWe thank the reviewers again for their thorough and constructive comments, which have\nhelped us assemble an improved manuscript.\nReferences\nAbrahams, B. S., D. E. Arking, D. B. Campbell, H. C. MeXord, E. M. Morrow, L. A. Weiss, I.\nMenashe, T. Wadkins, S. Banerjee-Basu and A. Packer (2013). \"SFARI Gene 2.0: A\ncommunity-driven knowledgebase for the autism spectrum disorders (ASDs).\" Molecular\nAutism 4(1): 1-3.\nBurgaya, F., M. Toutant, J. M. Studler, A. Costa, M. Le Bert, M. Gelman and J. A. Girault\n(1997). \"Alternatively spliced focal adhesion kinase in rat brain with increased\nautophosphorylation activity.\" J Biol Chem 272(45): 28720-28725.\nCharlet-Berguerand, N., H. Le Hir, M. Incoronato, U. di Porzio, Y. Yu, S. Jing, V. de Franciscis\nand C. Thermes (2004). \"Expression of GFRalpha1 receptor splicing variants with diXerent\nbiochemical properties is modulated during kidney development.\" Cell Signal 16(12): 1425-\n1434.\nDixit, A., O. Pamas, B. Y. Li, J. Chen, C. P. Fulco, L. Jerby-Amon, N. D. Marjanovic, D. Dionne,\nT. Burks, R. Raychowdhury, B. Adamson, T. M. Norman, E. S. Lander, J. S. Weissman, N.\nFriedman and A. Regev (2016). \"Perturb-Seq: Dissecting Molecular Circuits with Scalable\nSingle-Cell RNA Profiling of Pooled Genetic Screens.\" Cell 167(7): 1853-+.\nGantner, C. W., C. P. J. Hunt, J. C. Niclis, V. Penna, S. J. McDougall, L. H. Thompson and C. L.\nParish (2021). \"FGF-MAPK signaling regulates human deep-layer corticogenesis.\" Stem Cell\nReports 16(5): 1262-1275.\nGarcia-Cabau, C., A. Bartomeu, G. Tesei, K. C. Cheung, J. Pose-Utrilla, S. Picó, A.\nBalaceanu, B. Duran-Arqué, M. Fernández-Alfara, J. Martín, C. De Pace, L. Ruiz-Pérez, J.\nGarcía, G. Battaglia, J. J. Lucas, R. Hervás, K. LindorX-Larsen, R. Méndez and X. Salvatella\n(2024). \"Mis-splicing of a neuronal microexon promotes CPEB4 aggregation in ASD.\" Nature\n2024: 1-8.\nGonatopoulos-Pournatzis, T., M. Aregger, K. R. Brown, S. Farhangmehr, U. Braunschweig, H.\nN. Ward, K. C. H. Ha, A. Weiss, M. Billmann, T. Durbic, C. L. Myers, B. J. Blencowe and J.\nMoXat (2020). \"Genetic interaction mapping and exon-resolution functional genomics with\na hybrid Cas9–Cas12a platform.\" Nature Biotechnology 38(5): 638-648.\nHan, H., A. J. Best, U. Braunschweig, N. Mikolajewicz, J. D. Li, J. Roth, F. Chowdhury, F.\nMantica, S. Nabeel-Shah, G. Parada, K. R. Brown, D. O'Hanlon, J. Wei, Y. Yao, A. A. Zid, L. C.\nComsa, M. Jen, J. Wang, A. Datti, T. Gonatopoulos-Pournatzis, R. J. Weatheritt, J. F.\nGreenblatt, J. L. Wrana, M. Irimia, A. C. Gingras, J. MoXat and B. J. Blencowe (2022).\n\"Systematic exploration of dynamic splicing networks reveals conserved multistage\nregulators of neurogenesis.\" Mol Cell 82(16): 2982-2999 e2914.\nJiang, L. D., C. Dalgarno, E. Papalexi, I. Mascio, H. H. Wessels, H. Y. Yun, N. Iremadze, G.\nLithwick-Yanai, D. Lipson and R. Satija (2025). \"Systematic reconstruction of molecular\n24\npathway signatures using scalable single-cell perturbation screens.\" Nature Cell Biology\n27(3).\nLin, L., M. Zhang, P. Stoilov, L. Chen and S. Zheng (2020). \"Developmental Attenuation of\nNeuronal Apoptosis by Neural-Specific Splicing of Bak1 Microexon.\" Neuron 107(6): 1180-\n1196 e1188.\nMitra, S. K., D. A. Hanson and D. D. Schlaepfer (2005). \"Focal adhesion kinase: in command\nand control of cell motility.\" Nature Reviews Molecular Cell Biology 2005 6:1 6(1): 56-68.\nPapalexi, E., E. P. Mimitou, A. W. Butler, S. Foster, B. Bracken, W. M. Mauck, H. H. Wessels,\nY. H. Hao, B. Z. Yeung, P. Smibert and R. Satija (2021). \"Characterizing the molecular\nregulation of inhibitory immune checkpoints with multimodal single-cell screens.\" Nature\nGenetics 53(3): 322-+.\nParatcha, G., F. Ledda and C. F. Ibáñez (2003). \"The neural cell adhesion molecule NCAM is\nan alternative signaling receptor for GDNF family ligands.\" Cell 113(7): 867-879.\nParras, A., H. Anta, M. Santos-Galindo, V. Swarup, A. Elorza, J. L. Nieto-Gonzalez, S. Pico, I.\nH. Hernandez, J. I. Diaz-Hernandez, E. Belloc, A. Rodolosse, N. N. Parikshak, O.\nPenagarikano, R. Fernandez-Chacon, M. Irimia, P. Navarro, D. H. Geschwind, R. Mendez\nand J. J. Lucas (2018). \"Autism-like phenotype and risk gene mRNA deadenylation by\nCPEB4 mis-splicing.\" Nature 560(7719): 441-446.\nPatowary, A., P. Zhang, C. Jops, C. K. Vuong, X. Z. Ge, K. C. Hou, M. Kim, N. H. Gong, M.\nMargolis, D. Vo, X. S. Wang, C. Y. Liu, B. Pasaniuc, J. J. Li, M. J. Gandal and L. de la Torre-\nUbieta (2024). \"Developmental isoform diversity in the human neocortex informs\nneuropsychiatric risk mechanisms.\" Science 384(6698).\nSantinha, A. J., E. Klingler, M. Kuhn, R. Farouni, S. Lagler, G. Kalamakis, U. Lischetti, D.\nJabaudon and R. J. Platt (2023). \"Transcriptional linkage analysis with in vivo AAV-Perturb-\nseq.\" Nature 2023 622:7982 622(7982): 367-375.\nSatterstrom, F. K., J. A. Kosmicki, J. Wang, M. S. Breen, S. De Rubeis, J. Y. An, M. Peng, R.\nCollins, J. Grove, L. Klei, C. Stevens, J. Reichert, M. S. Mulhern, M. Artomov, S. Gerges, B.\nSheppard, X. Xu, A. Bhaduri, U. Norman, H. Brand, G. Schwartz, R. Nguyen, E. E. Guerrero,\nC. Dias, C. Autism Sequencing, P.-B. C. i, C. Betancur, E. H. Cook, L. Gallagher, M. Gill, J. S.\nSutcliXe, A. Thurm, M. E. Zwick, A. D. Borglum, M. W. State, A. E. Cicek, M. E. Talkowski, D.\nJ. Cutler, B. Devlin, S. J. Sanders, K. Roeder, M. J. Daly and J. D. Buxbaum (2020). \"Large-\nScale Exome Sequencing Study Implicates Both Developmental and Functional Changes in\nthe Neurobiology of Autism.\" Cell 180(3): 568-584 e523.\nToutant, M., A. Costa, J. M. Studler, G. Kadare, M. Carnaud and J. A. Girault (2002).\n\"Alternative splicing controls the mechanisms of FAK autophosphorylation.\" Mol Cell Biol\n22(22): 7731-7743.\nWang, J., Y. Jia, S. Zhao, X. Zhang, X. Wang, X. Han, Y. Wang, M. Ma, J. Shi and L. Liu (2017).\n\"BIN1 reverses PD-L1-mediated immune escape by inactivating the c-MYC and\nEGFR/MAPK signaling pathways in non-small cell lung cancer.\" Oncogene 2017 36:45\n36(45): 6235-6243.\nWang, W., L. Huang, Y. Huang, J. w. Yin, A. J. Berk, J. M. Friedman and G. Wang (2009).\n\"Mediator MED23 links insulin signaling to the adipogenesis transcription cascade.\"\nDevelopmental cell 16(5): 764-771.\nWang, W. J., J. C. Kuo, W. Ku, Y. R. Lee, F. C. Lin, Y. L. Chang, Y. M. Lin, C. H. Chen, Y. P.\nHuang, M. J. Chiang, S. W. Yeh, P. R. Wu, C. H. Shen, C. T. Wu and R. H. Chen (2007). \"The\ntumor suppressor DAPK is reciprocally regulated by tyrosine kinase Src and phosphatase\nLAR.\" Molecular cell 27(5): 701-716.\nWang, X., L. P. Weng and Q. Yu (2000). \"Specific inhibition of FGF-induced MAPK activation\nby the receptor-like protein tyrosine phosphatase LAR.\" Oncogene 2000 19:19 19(19): 2346-\n2353.\nWood, S. N. (2011). \"Fast stable restricted maximum likelihood and marginal likelihood\nestimation of semiparametric generalized linear models.\" Journal of the Royal Statistical\nSociety Series B-Statistical Methodology 73: 3-36.\nYang, X., M. Zhao, M. Xia, Y. Liu, J. Yan, H. Ji and G. Wang (2012). \"Selective requirement for\nMediator MED23 in Ras-active lung cancer.\" Proceedings of the National Academy of\nSciences of the United States of America 109(41).\nYang, Y., C. Li, Z. Chen, Y. Zhang, Q. Tian, M. Sun, S. Zhang, M. Yu and G. Wang (2023). \"An\nintellectual disability-related MED23 mutation dysregulates gene expression by altering\nchromatin conformation and enhancer activities.\" Nucleic acids research 51(5): 2137-\n2150.\n26"}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this revision, the authors have clarified some challenging concepts. For instance, adoption of differentiation rank and eCDFs is more reader-friendly and illustrates the nature of gene expression changes upon microexon knockouts, i.e., precocious expression of neural genes. They have also addressed some points about the deletion of constitutive exons from the same set of host genes. Below are additional comments on the revised manuscript the authors might consider:\n1. On microexon deletion efficiency in Fig. 1: We do not fully agree with the authors that the emphasis of Fig. 1 is to optimize the CHyMErA system by comparing editing efficiencies between Cas12a constructs. The representative efficiency of the system is a reliability indicator for all downstream efforts. There is still some concern regarding whether the deletion efficiencies listed in Fig. 1c-d reflect clonal ground truth. In Revision Fig. 1b, % edited clones (one assumes) provides a binary measure of deletions (heterozygous and homozygous deletions), and does not provide a clear indication of how frequent homozygous deletions may be, or account for the wider set of possibilities in bulk populations (het. versus hom. deletions). For example, 5/5 clonal lines with heterozygous deletions may be described as 100 % edited clones, but in the bulk population the same clones will have 50% deletions. Moreover, as the length normalized band intensities only corrects for differences in intercalations, it does not correct for the more concerning amplification bias of shorter deletion amplicons. It seems disingenuous to ignore shorter amplicon bias. Perhaps linkage is lost in ddPCRs, but it could still provide a more grounded estimate of bulk editing. Beyond further efforts to resolve true editing efficiency, the authors may consider listing this as a limitation.\n2. The notion that heterozygote microexon deletions could have phenotypes is plausible, but in this case the clonal het. lines should be used to demonstrate transcriptomic effects for validation targets.\n3. The reliability of QDA still remains to be determined. The authors point to validation of example target exons in this manuscript, but this is unsatisfactory. New analytic pipelines such as this should be benchmarked against published datasets. A good example is single cell CRISPR knock out experiments where it is possible to quantify depletion of target gene.\n4. How many KO cells/cluster or subcluster per target for gene expression analyses?\n5. The representation of guide pairs in the sequenced nuclei library (and potential sources of bias) still needs further illustration. While we understand some basic ideas about #cells/target a detailed description of what the identities of these cells are will be most informative to the reader. In particular, as this study performs DGE on KO v non-targeting subclusters in pseudo-bulk, it is essential to understand how many cells/target microexon KO within each subcluster. We have noted – with concern – that there is only a moderate correlation between bulk and snRNA comparisons (including Figs. 2g v 3c Gfra1 KO and Ext Fig. 5c ) and questioned the DGE outcome of the screen. Understanding whether low Ns could be a contributing weakness is quite important.\n6. The microexon KO cells are broadly located in subclusters 0,2,4 and 11, all of which represent neuroepithelial (NE) identities (Fig. 2d). There is a strong, general trend for microexon KO cells to be found in neuroepithelial state, with most non-targeted cells in neuron subclusters (orange points in Fig. 2c). The authors have clarified that the microexons are included in NE stage, but it is not clear if KO cell enrichment in NE has true biological meaning, i.e., differentiation defects, or if this is some unexpected technical issue. As there is no replicate screen or characterization of differentiation using clonal microexon KO cells, the robustness of these findings remains uncertain.\n7. Related to Fig. 2g, it would be useful to understand which subclusters enrich precocious neurogenesis signatures.\n8. It is true that “estimating the ‘true’ proportion of edited/effected cells in bulk cell populations represents a more general challenge encountered in CRISPR-based screens”, however there are general approaches to mitigate concerns, for example independent replicate screens which have not been performed.\n9. On modest correlation between changes in gene expression between snRNA and bulk RNAseq datasets: technical differences between the modalities should be described using comparisons with data generated in this study. One cannot simply take the authors word that there are big differences modalities. Extended Fig. 5c shows an alarmingly poor correlation between bulk RNA-seq from a clonal line and CHyMErA-seq. GO-term overlaps as shown in Fig. 3b could indicate bona fide biology or easy-to-obtain terms within the context of neuronal differentiation. As this one of few validations done in this study, it highlights some critical challenges that remain with demonstrating robustness of screening with CHyMErA system.\n10. No experiments were done to further validate transcriptomic or other phenotypes in clonal KO lines.\n(Remarks on code availability)", "reviewer_2": "The authors have addressed all the points I had raised in the original review and I have no further comments.\n(Remarks on code availability)"}, "author_response": "Reviewer #1 (Remarks to the Author):\nIn this revision, the authors have clarified some challenging concepts. For instance, adoption of di9erentiation rank and eCDFs is more reader-friendly and illustrates the nature of gene expression changes upon microexon knockouts, i.e., precocious expression of neural genes. They have also addressed some points about the deletion of constitutive exons from the same set of host genes. Below are additional comments on the revised manuscript the authors might consider:\nWe thank the reviewer for their additional helpful comments on our manuscript. Our\nresponses to these are provided below, along with a summary of the corresponding changes\nwe have made to our manuscript (with all new text changes highlighted in yellow).\n1. On microexon deletion e9iciency in Fig. 1: We do not fully agree with the authors that the\nemphasis of Fig. 1 is to optimize the CHyMErA system by comparing editing e9iciencies\nbetween Cas12a constructs. The representative e9iciency of the system is a reliability\nindicator for all downstream e9orts. There is still some concern regarding whether the\ndeletion e9iciencies listed in Fig. 1c-d reflect clonal ground truth. In Revision Fig. 1b, %\nedited clones (one assumes) provides a binary measure of deletions (heterozygous and\nhomozygous deletions), and does not provide a clear indication of how frequent\nhomozygous deletions may be, or account for the wider set of possibilities in bulk\npopulations (het. versus hom. deletions). For example, 5/5 clonal lines with heterozygous\ndeletions may be described as 100 % edited clones, but in the bulk population the same\nclones will have 50% deletions.\nWe apologize for the confusion. In our previous Rebuttal Fig. 1b, we were in fact showing the\npercent of edited alleles across clonal cell lines. Therefore, ‘% Edited Clones’ indicates the\npercentage of editing while considering the zygosity of the deletion; i.e. 5/5 homozygous\ndeletions would be considered 100% edited, whereas 5/5 heterozygous deletions would be\nconsidered 50% edited.\nMoreover, as the length normalized band intensities only corrects for di9erences in\nintercalations, it does not correct for the more concerning amplification bias of shorter\ndeletion amplicons. It seems disingenuous to ignore shorter amplicon bias. Perhaps linkage\nis lost in ddPCRs, but it could still provide a more grounded estimate of bulk editing. Beyond\nfurther e9orts to resolve true editing e9iciency, the authors may consider listing this as a\nlimitation.\nWe agree with the reviewer that there could be potential amplification bias favouring\ndetection of edited versus unedited bands. However, as previously shown in our original\nrebuttal, the editing e9iciencies we observe are consistent with those observed in clonal\npopulations, which serve as a ‘ground-truth’ for estimating editing e9iciency.\nWe appreciate the reviewer’s repeated suggestion to apply ddPCR to address this issue.\nHowever, in addition to our previous response regarding the inability of this approach to\naddress the stated issues (i.e. zygosity of deletions, which we nevertheless have addressed\nby analysing clonal lines), we would also like to point out that this approach is far from\nstraightforward and would require the design and optimization of distinct sets of primers\nwith comparable e9iciencies to capture edited and unedited genomes.\nIn order to further address the reviewer’s concern, we have now provided an independent\nmeasure of the editing e9iciency of our enhanced CHyMErA system using Oxford nanopore\nsequencing of amplified genomic regions targeted for deletion. This approach allows an\nassessment of editing e9iciency independently of possible length-dependent measurement\nbiases associated with agarose gel densitometry, such as ethidium bromide intercalation.\nWe observe significantly reduced sequencing coverage over the expected deletion region for\nanalyzed constructs, relative to the non-targeting control (Rebuttal Fig. 1a). We calculated\nthe editing e9iciency in each sample based on reads spanning Cas9/Cas12a deletion sites\nand the intervening region and observed 45% and 52% editing e9iciency for constructs 1 and\n2, respectively (Rebuttal Fig. 1b), whereas only 2% of non-targeting reads were classified as\n‘edited’ (although likely due to misclassified, truncated reads). The values for editing\ne9iciency determined using Nanopore sequencing are comparable (i.e. within 10%) to\nestimates of editing e9iciencies calculated using semi-quantitative measurements of\nagarose gel band intensities. This validation approach is now mentioned in the Methods.\nPage 21:\n“To independently determine approximate editing e2iciencies of the eCHyMErA method,\namplified PCR products shown in Figure 1c were analyzed by Oxford Nanopore long-read\nsequencing (performed by Plasmidsaurus Inc.). Raw reads were aligned to sequences\ncorresponding to the amplified regions. Editing e2iciency was quantified by comparing the\nrelative proportions of reads mapping to Cas9/Cas12a deleted regions versus adjacent\nunedited sites. Editing e2iciencies estimated by quantification of gel band intensities and\nlong-read sequencing were within 10% of each other.”\nRebuttal Figure 1. Quantifying CHyMErA-directed deletions by Oxford Nanopore long-read sequencing. a)\nRead coverage surrounding the Ptk2 microexon ‘A’ targeted for CHyMErA deletion. Cut sites for Cas9 and\nCas12a are indicated by red and blue dashed lines, respectively; position of microexon is indicated by\nyellow solid lines. b) Quantification of reads representing unedited (grey), repair (orange) or deletion (red)\nof the intervening genomic regions. Overall read counts analyzed: non-targeting (n=719 reads), Ptk2\ndeletion construct 1 (n = 624 reads) and Ptk2 deletion construct 2 (n=344 reads).\nIn summary, by determining ground truth editing e9iciencies through the time-consuming\ngeneration and genotyping of clonal cell lines, and by using two independent methods for\nquantification of CHyMErA editing e9iciencies in these bulk cell populations, we observe\nsimilar e9iciencies. Furthermore, it is worth noting that we and others have extensively\ncharacterized and validated the CHyMErA system in previous publications directed at\nmeasuring phenotypes in response to deletions of targeted exons1 , 2. This includes a recent\npublication2 that validates CHyMErA-detected exon deletion growth phenotypes using base-\nediting-induced skipping. To document how we have quantified editing e9iciency using two\nindependent types of measurements, and to acknowledge limitations in the accuracy of\nestimating editing e9iciencies, we have revised the Figure 1 legend:\nFigure 1c.\n“Editing e2iciencies were estimated from length-normalized quantification of ratios of edited\nand total band intensities, and are consistent (i.e. within 10%) of estimates determined by\nquantification of Oxford Nanopore sequencing reads mapping specifically to unedited and\nedited genomic fragments (Methods).”\n2. The notion that heterozygote microexon deletions could have phenotypes is plausible, but\nin this case the clonal het. lines should be used to demonstrate transcriptomic e9ects for\nvalidation targets.\nWe have already addressed this point by including validation experiments in our manuscript\nin which dCasRx-RBM25 was used to directed modulation of microexon splicing.\nSpecifically, targeting the Clasp1(A) microexon results in partial (~50%) exon skipping, yet\nthis level is su9icient to recapitulate transcriptomic phenotypes observed in the screen. As\nalready mentioned in our response, previous focused analyses comparing heterozygous and\nhomozygous deletions of specific neural-di9erential exons have similarly resulted in\nsignificant yet reduced transcriptomic and additional phenotypes for the heterozygous\ndeletions3, 4. We have clarified this point in the Discussion.\nPage 16:\n“Moreover, consistent with previous studies demonstrating that heterozygous deletion of\nalternative exons can have a significant impact on transcriptomic regulation3, 4, the results of\nour CHyMErA-seq screen and follow-up experiments using dCasRx-RBM25-directed exon\nmanipulation further demonstrate that induced partial skipping of microexons can lead to\nsignificant phenotypes, an observation that may also have relevance for interpreting the\nsignificance of partial exon skipping levels detected in patient brain samples.”\n3. The reliability of QDA still remains to be determined. The authors point to validation of\nexample target exons in this manuscript, but this is unsatisfactory. New analytic pipelines\nsuch as this should be benchmarked against published datasets. A good example is single\ncell CRISPR knock out experiments where it is possible to quantify depletion of target gene.\nWe thank the reviewer for this excellent suggestion. We have performed the requested\nbenchmarking using a previously published single-cell perturbation screen that utilized\nCas9-mediated knockout of target genes5. As requested, we scored all cells according to\ntheir QDA-derived perturbation signature and quantified di9erences in the expression level\nof the target gene in cells classified as ‘perturbed’, relative to control cells expressing a non-\ntargeting guide. We observe that 12/15 perturbations for which we detect a transcriptomic\nsignature also have significant reductions in expression levels of the targeted gene (P < 0.05,\nBonferroni-corrected Wilcoxon rank-sum test) (Rebuttal Fig. 2). This benchmarking\ncomparison is now mentioned on page 7:\n“Supporting the reliability of these results, benchmarking of our analysis pipeline against a\npublished single-cell CRISPR screen dataset5 reveals that 12/15 (80%) of perturbations\ncalled by our method have significant levels of reduction in expression of the target gene (P\n< 0.05, Bonferroni-corrected Wilcoxon rank-sum test).”\nRebuttal Figure 2. Expression of target gene for non-targeting cells (grey) and cells classified as perturbed\nby QDA pipeline (red). Data from Papalexi et al 2021. *P < 0.05, **P < 0.01, ***P < 0.001. Bonferroni-\ncorrected Wilcoxon rank-sum test.\n4. How many KO cells/cluster or subcluster per target for gene expression analyses?\nThis information is provided in the legend for Figure 2d:\n“d) Enrichment of ‘knockout’ cells compared to non-targeting cells within each subcluster,\ncorresponding to deletion of Ralgapb (n = 105), Ptprf (microexon A; n = 158), Ptk2\n(microexons A, B and C; n = 112), Med23 (n = 91), Gfra1 (n = 103), Clasp1 (microexon A; n =\n130), Camta1 (n = 34) and Bin1 (n = 180) microexons (bottom).”\n5. The representation of guide pairs in the sequenced nuclei library (and potential sources\nof bias) still needs further illustration. While we understand some basic ideas about\n#cells/target a detailed description of what the identities of these cells are will be most\ninformative to the reader. In particular, as this study performs DGE on KO v non-targeting\nsubclusters in pseudo-bulk, it is essential to understand how many cells/target microexon\nKO within each subcluster.\nTo address this concern, we have plotted the number of cells represented by each\nmicroexon-targeting guides for each cluster and added this to Extended Data Figure 3c. This\nnew plot is referenced on page 7 of the revised text.\n“We clustered cells and analysed marker gene expression within these trajectories to reveal\nthe presence of expected neuroepithelial, radial glial, and neuronal clusters within the\nneuronal lineage (Fig. 2b, Extended Data Figs. 3a and 3b), with cells expressing microexon\ntargeting guides for 24/36 detected microexon targets represented by at least 25 cells in each\ncluster within the neuronal lineage (Extended Data Fig. 3c).”\nRebuttal Figure 3. Number of cells assigned to each microexon target according to cluster.\nWe have noted – with concern – that there is only a moderate correlation between bulk\nand snRNA comparisons (including Figs. 2g v 3c Gfra1 KO and Ext Fig. 5c) and questioned\nthe DGE outcome of the screen. Understanding whether low Ns could be a contributing\nweakness is quite important.\nWe acknowledge the reviewer’s repeated concerns regarding the low correlation between\nthe single-nuclei and bulk RNA-seq signatures. However, we would like to again emphasize\nthat there are major inherent di9erences between these methods that will confound deriving\ndirect correlations. However, we have observed and reported positive and statistically\nsignificant correlations between the two methods, despite the manifold possible\nconfounding factors.\nTo further illustrate the challenges, we demonstrate below the substantial di9erences\nbetween the transcriptomic data collected from bulk, whole-cell RNA-seq and single-\nnuclear RNA seq, when calculating correlations between gene expression patterns from\nthese data. As expected, there are no OVERALL significant correlations between expression\nlevels measured using single-nuclear RNA-seq and the bulk RNA-seq, yet, despite this, we\nstill observe positive and significant correlations for microexon perturbation-dependent\ngene expression signatures. (Rebuttal Fig. 4).\nsample_type\nNT 6\nNT 1\nNT 3\nNT 2\nNT 7 Correlation sample_type\nNT 8 1 Bulk RNA−seq\n0.5 snRNA−seq\nNT 4\nNT DIV7 B3 R2 −1\nNT DIV7 B4 R1\nNT DIV7 B4 R2\nNT ESC B3 R1\nNT ESC B4 R2\nNT\nR1\nB3\nR2\nB4\nR1\nB4\nR2\nB3\nR1\nB4\nR2\nNT\nDI V\nN\nDI V\nN\nDI V\nN\nDI V\nN\nE S C\nN T\nE S C\nRebuttal Figure 4. Visualizing the correlation between gene expression for non-targeting cells in single-\nnuclear RNA-seq (pink) and whole cell bulk RNA-seq (blue).\n6. The microexon KO cells are broadly located in subclusters 0,2,4 and 11, all of which\nrepresent neuroepithelial (NE) identities (Fig. 2d). There is a strong, general trend for\nmicroexon KO cells to be found in neuroepithelial state, with most non-targeted cells in\nneuron subclusters (orange points in Fig. 2c). The authors have clarified that the microexons\nare included in NE stage, but it is not clear if KO cell enrichment in NE has true biological\nmeaning, i.e., di9erentiation defects, or if this is some unexpected technical issue. As there\nis no replicate screen or characterization of di9erentiation using clonal microexon KO cells,\nthe robustness of these findings remains uncertain.\nContrary to the reviewer’s assertion, most of the non-targeting cells are in fact not located in\nthe neuronal subclusters. Rather, they are dispersed throughout the clusters (see Rebuttal\nFigure 5).\nRebuttal Figure 5. Bar plot of relative enrichment of non-targeting cells in each of the clusters relative to\nthe overall proportion of cells. Error bars indicate 95% Wilson score interval.\nWe are puzzled by the reviewer’s suggestion that we have not further characterized this\nphenotype, since we have indeed characterized clonal cell lines for Gfra1 microexon\ndeletion, as well phenotypes using multiple dCasRx-RBM25 cell lines designed to validate\nthe e9ects of induced microexon skipping on di9erentiating cells (please see Figures 3 and\n4 and the corresponding text).\n7. Related to Fig. 2g, it would be useful to understand which subclusters enrich precocious\nneurogenesis signatures.\nFigure 2g illustrates the distribution of di9erentially expressed genes resulting from\nmicroexon deletion with respect to the neuron di9erentiation score. As previously described,\nthese scores represent the correlation between gene expression and neuron lineage\nprogression. The information requested by the reviewer is in Figure 2d, where we show which\nsubclusters are enriched in cells classified as ‘perturbed’.\n8. It is true that “estimating the ‘true’ proportion of edited/e9ected cells in bulk cell\npopulations represents a more general challenge encountered in CRISPR-based screens”,\nhowever there are general approaches to mitigate concerns, for example independent\nreplicate screens which have not been performed.\nIt is unclear to us how a replicate screen would elucidate the ‘true’ portion of edited cells.\nFurthermore, as described in the previous rebuttal, we are performing extensive internal\nreplication by applying nine distinct guide pairs to target the same microexons. This strategy\nis significantly more cost e9icient than performing multiple replicate screens with fewer\nguides per target (which would be required to maintain coverage while sequencing the same\ntotal number of cells) and avoids encountering batch e9ects that might otherwise\ncompromise making direct comparisons between replicates.\n9. On modest correlation between changes in gene expression between snRNA and bulk\nRNAseq datasets: technical di9erences between the modalities should be described using\ncomparisons with data generated in this study. One cannot simply take the authors word\nthat there are big di9erences modalities.\nPlease refer above to rebuttal Figure 4 and the associated description of this analysis.\nExtended Fig. 5c shows an alarmingly poor correlation between bulk RNA-seq from a clonal\nline and CHyMErA-seq. GO-term overlaps as shown in Fig. 3b could indicate bona fide\nbiology or easy-to-obtain terms within the context of neuronal di9erentiation. As this one of\nfew validations done in this study, it highlights some critical challenges that remain with\ndemonstrating robustness of screening with CHyMErA system.\nWe hope that by illustrating the major di9erences between these approaches (Rebuttal Fig.\n4) the reviewer can appreciate the challenges when making direct comparisons between\nsingle-nuclear RNA-seq and whole cell bulk RNA-seq. However, as pointed out above and\non pages 10 and 11 of our manuscript, despite these di9erences, positive and significant\ncorrelations were observed between altered gene expression signatures detected using\nthese methods, in addition to the significant overlap in GO terms. Together we believe these\nresults strongly support the validity of the phenotypic e9ects we observe in the CHyMErA-\nseq screen. Furthermore, the GO terms are all scored against a background of genes that are\nexpressed in our in vitro di9erentiated neurons and therefore do not represent ‘easy-to-\nobtain’ enrichments. Even so, we followed-up on these GO-term enrichments, and in the\ncase of Gfra1, when using clonal lines, we observe that in addition to displaying an\nenrichment of gene expression signatures that correspond to pathways involved in cell\nmorphogenesis and neuron projection formation (Fig. 3b), we observe that that these\nneurons display morphological features that reflect the transcriptional changes (i.e. longer\nneurites) (Figs. 3d and 3e).\n10. No experiments were done to further validate transcriptomic or other phenotypes in\nclonal KO lines.\nThis statement is incorrect. As clarified above and in our manuscript, we have validated\nsome or our observations using clonal lines in which microexons are targeted at the RNA-\nlevel using dCasRx-RBM25 (please refer to Figures 3 and 4 and pages 9-12).\nReferences\n1. Gonatopoulos-Pournatzis, T. et al. Genetic interaction mapping and exon-resolution\nfunctional genomics with a hybrid Cas9-Cas12a platform. Nat Biotechnol 38, 638-\n648 (2020).\n2. Xiao, M.S. et al. Genome-scale exon perturbation screens uncover exons critical for\ncell fitness. Mol Cell (2024).\n3. Gueroussov, S. et al. An alternative splicing event amplifies evolutionary di9erences\nbetween vertebrates. Science 349, 868-873 (2015).\n4. Han, H. et al. Systematic exploration of dynamic splicing networks reveals conserved\nmultistage regulators of neurogenesis. Mol Cell 82, 2982-2999 e2914 (2022).\n5. Papalexi, E. et al. Characterizing the molecular regulation of inhibitory immune\ncheckpoints with multimodal single-cell screens. Nat Genet 53, 322-331 (2021)."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I appreciate the authors' efforts to clarify my previous comments. I have no more additional points.\n(Remarks on code availability)\nN/A"}, "author_response": ""}]

Decision letter