Animals
All animals used in this study were approved by the Institutional Animal Care and Use Committee at ISTA and by the Bundesministerium für Bildung, Wissenschaft und Forschung, Austria (approval numbers BMWFW−66.018/0008-WF/II/3b/2014, BMWFW-66.018/0012-WF/V/3b/2015, BMWFW-66.018/0012-WF/V/3b/2017, BMWFW-66.018/0015-V/3b/2019, BMWFW-66.018/0032-V/3b/2019, BMBWF-V/3b/2020-0.342.159, BMBWF-V/3b/2020-0.148.791, BMBWF-V/3b/2021-0.291.172, BMBWF-V/3b/2021-0.291.177, BMBWF-V/3b/2022-0.292.788 and BMBWF-V/3b/2022-0.121.445). All our mutant lines were maintained in a C57BL/6J background. To maintain this background, colonies are continuously refreshed by mating mutants with C57BL/6J animals drawn from a continuously refreshed central C57BL/6J colony. All experiments were performed on mice ranging from E14.5 to P4 and P14. Embryonic time points were determined by plug checks after timed matings, defining E0.5 as the morning post-coitum. Animals were kept in the Preclinical Facility at ISTA, housed in commercially available individually ventilated cages under defined standard laboratory conditions (room temperature 22 ± 1 °C, relative humidity 55 ± 10%) on a 12-h light–dark cycle (lights on at 07:00). Animals were housed in groups of 3–4 animals per cage, with food and water available ad libitum. Experiments were carried out under specific pathogen-free conditions, and the health status of the mouse lines was routinely checked by a veterinarian. All transgenic mouse lines were backcrossed into a C57BL/6J background a minimum of two times. Both females and males were used for experiments.
Generation of constitutive heterozygous mouse lines
Constitutive heterozygous mice (that is, Ash1l+/−, Kdm6b+/−, Kmt5b+/−, Trip12+/−, Usp7+/− and Wac+/−) were generated within the Transgenic Unit at IST Austria in a C57BL/6J background by zygote electroporation of the Cas9 protein (634621, TakaraBio) and custom-made sgRNAs. To generate Ash1l+/− mice, one sgRNA (CAATCACCATTCCGCTATCC AGG) targeting exon 12 of the Ash1l gene was used. For Kdm6b+/−, two sgRNAs (CCCAAGGACCCGAGGTGATA GGG and TACGTATGAGGAGCGAACCC TGG) targeting intronic regions flanking exon 19 of the Kdm6b gene were used. For Kmt5b+/−, two sgRNAs (CCCTCACGACCCTACTACTG TGG and TTTACGTCTCAAGTCACACT GGG) targeting intronic regions flanking exon 9 of the Kmt5b gene were used. Similarly, to generate Trip12+/− mice, two sgRNAs (AGCGTATTCTCACTTTGATC TGG and TAAATTTGTAGAGACAGTTA AGG) targeting intronic regions flanking exon 3 of the Trip12 gene were used. Similarly, for Usp7+/−, two sgRNAs (CTCTCCCGAACATTATGAAG AGG and TGCCACGTCTTATACACATT AGG) targeting intronic regions flanking exon 6 of the Usp7 gene were used. To generate Wac+/− mice, one sgRNA (ATGTCGGAGGCGCAGATGTA GGG) targeting exon 8 of the Wac gene was used.
Superovulation was performed to produce large numbers of zygotes. In brief, 3–4-week-old donor females (in-house) were hormone primed by intraperitoneal injection of 5 IU of pregnant mare serum gonadotropin. Of human chorionic gonadotropin, 5 IU was administered at 46–48 h after pregnant mare serum gonadotropin treatment, followed by 1:1 mating with stud males (in-house). Zygotes were isolated from the oviduct of mated females at 0.5 days post-coitum. M2 with hyaluronidase (MR-051-F, Sigma) was used to dissociate cumulus cells. Zygotes were then cultured in small drops of KSOM (MR-121-D, Sigma) under mineral oil (M-5310, Sigma) at 37 °C and 5% CO2, until the electroporation. Electroporation of mouse embryos was performed similarly as previously described63 with the adaptation of using Opti-MEM (31985062, Thermo Fisher) instead of PBS. Finally, F0 founder mice were genotyped using standard PCR amplification with primers 500 bp to 1 kb away from guide RNA-binding sites. F1 mice were generated by backcrossing founder mice with C57BL/6J animals. The F1 mice were again genotyped by PCR, and positive amplicons covering the deleted sequence were Sanger sequenced to confirm precise targeting. All transgenic mouse lines were back-crossed into a C57BL/6J background a minimum of two times, and predicted off-targets were checked. A list of oligonucleotides used is provided in Supplementary Table 1. The constitutive knockout mouse line (HnrnpU+/−) was obtained by mating HnrnpUWT/flox (Hnrnpu<tm1.1Tman>/J, strain: #032187, RRID:IMSR_JAX:032187) mice with a CMV-Cre line (B6.C-Tg(CMV-cre)1Cgn/J; strain #:006054, RRID:IMSR_JAX:006054) and subsequent backcrossings to C57BL/6J WT animals. Crossing of HnrnpUWT/flox mice with homozygous animals of the CMV-CreCre/Cre line results in ubiquitous deletion of one allele of HnrnpU spanning exons 4–14.
Single-nucleus multi-omics profiling of gene expression (snRNA-seq) and chromatin accessibility (snATAC-seq) from the same nucleus (snMultiome-seq)
Sample preparation and nuclei isolation
All the sequenced mutant animals were generated by mating a C57BL/6J female from the continuously refreshed central C57BL/6J colony with a mutant male maintained on the same C57BL/6J background (for example, C57BL/6J female × Ash1l+/− male = litter sequenced). The only exceptions were homozygous Bckdk mice whose heterozygous breeders were regularly refreshed with the same central C57BL/6J colony. Control animals were C57BL/6J × C57BL/6J offspring sourced from the same central colony. Dams used for mutant and control litters were drawn from the same central-colony cohorts and often littermates; at minimum, they were not separated by generations of inbreeding.
Brain tissue (forebrain and cerebellum) from E14.5, P4 and P14 mice was used for snMultiome-seq experiments. Animals coming from at least three separate litters were collected to avoid potential litter-specific biases. For the embryonic time point (that is, E14.5; n = 3 per genotype + sex), timed-pregnant females were decapitated, and embryos were rapidly dissected on ice. Cortical tissue was dissected in ice-cold PBS, meninges were removed, and forebrain structures were snap-frozen in liquid nitrogen and stored at −80 °C until nuclei isolation.
For postnatal time points (that is, P4 and P14; n = 3 per genotype + sex + time point), mice were decapitated, the brain was rapidly dissected on ice and subcortical structures were removed. The remaining forebrain (cortex + hippocampus) and cerebellum were separately snap-frozen in liquid nitrogen and stored at −80 °C until nuclei isolation.
Frozen tissue samples were dissociated to generate single-nuclei gene expression libraries following a modified protocol from 10X Genomics (Chromium Nuclei Isolation Kit CG000505, Rev A). Buffers were supplemented with RNase inhibitors (provided by 10X Genomics and 3335402001, Sigma). In brief, 200 µl of lysis buffer was added to the sample and dissociated with pestles provided by the vendor. After adding additional 300 µl of lysis buffer, the tissue was homogenized by pipette mixing and incubated on ice for 10 min. For the P14 forebrain, an additional 500 µl of EZ lysis buffer (NUC101, Sigma) was added to the samples to prevent clogging of the nuclei isolation columns. Afterwards, 500 µl of the cell suspension was loaded onto nuclei isolation columns and centrifuged for 20 s at 16,000 rcf at 4 °C, followed by vortexing the flow-through for 10 s and centrifugation for 3 min at 500 rcf at 4 °C in a swinging bucket rotor. The supernatant was removed, and the pellet was gently resuspended in 500 µl of debris removal buffer before centrifugation for 10 min at 700 rcf at 4 °C. After removal of the supernatant, the pellet was washed with 1 ml of wash buffer followed by centrifugation for 5 min at 500 rcf at 4 °C. For multiplexing, six samples were pooled after barcoding with CMOs (sequences of all CMOs used are provided in Supplementary Table 1). Each sample pellet was resuspended in a total of 180 µl of wash buffer containing 0.2 µM unique CMO anchor–barcode, vortexed for 3 s and incubated on ice for 5 min. CMO co-anchor was added to all samples of a multiplexing group with a final concentration of 0.2 µM, vortexed for 3 s and incubated on ice for 5 min. After adding 1 ml of wash buffer, samples were centrifuged for 5 min at 500 rcf at 4 °C. Depending on the tissue size, one (embryonic, P4 cortex and cerebellum) or two (P14 cortex) additional washing steps were performed. After removing the remaining supernatant, the pellet was resuspended in 50 µl of resuspension buffer. Nuclei were counted by flow cytometry (BD LSRFortessa) with Precision Count Beads (424902, BioLegend), and samples were diluted to obtain approximately 4,000 nuclei per sample, resulting in a total of 20,000 nuclei per lane of a 10X Chromium Controller (PN-1000204). GEMs were generated using Chromium Next GEM Chip J Single Cell kit (PN-1000230), Chromium Next GEM Single Cell Multiome GEM Kit A (PN-1000232) and Chromium Next GEM Single Cell Multiome Gel Bead Kit A (PN-1000231) as instructed by the manufacturer.
Library preparation and snMultiome sequencing
Library preparation was carried out according to the manufacturer’s protocol (Chromium Next GEM single cell multiome ATAC + gene expression, CG000338, Rev F) with the following modifications regarding the multiplexing reactions: during cDNA amplification, 0.2 µM ADT additive primer (5′-CCTTGGCACCCGAGAATTCC-3′) was added to the reaction. During cDNA cleanup, 120 µl of supernatant was saved to purify Hash-DNA using a 2.0× ratio of SPRIselect bead volume (B23319, Beckman Coulter) to sample volume. Hashtag libraries were generated using 5 ng of Hash-DNA, 1X KAPA2G Fast HotStart Readymix (7958935001, Roche) and unique combinations of 0.4 µM i5 and i7 indices (sequences are provided in Supplementary Table 1). Following amplification (98 °C for 2 min, 10 times (98 °C for 20 s, 54 °C for 30 s, 72 °C for 20 s), 72 °C for 5 min, 4 °C hold), libraries were cleaned up using a 1.2× ratio of SPRIselect bead volume (B23319, Beckman Coulter) to sample volume. Library quality and size were assessed with the High Sensitivity DNA Analysis Kit (5067-4626, Agilent Technologies) and the Bioanalyzer 2100. Libraries were sequenced by the Biomedical Sequencing Facility at the CeMM Research Center for Molecular Medicine of the Austrian Academy of Sciences, using the Illumina NovaSeq 6000 platform. For gene expression, libraries were pooled in a ratio of 83% RNA:17% HTO and sequenced on Illumina NovaSeq6000 systems using the sequencing format: 28 cycles for read 1, 10 cycles for i7, 10 cycles for i5 and 90 cycles for read 2. For snATAC-seq, libraries were sequenced on Illumina NovaSeq6000 systems using the sequencing format: 50 cycles for read 1N; 8 cycles for i7, 24 cycles for i5 and 49 cycles for read 2N.
Single-nucleus muli-omics data analysis
Pre-processing and demultiplexing, initial analysis and clustering
Samples used for multi-omics analyses were demultiplexed, and the raw sequencing data were converted to FASTQ format using bcl2fastq (v2.20.0.422). Joint processing of the gene expression and ATAC modalities was performed using the cellranger-arc software (v2.0.0). The mm10 reference genome was used for the alignment (refdata-cellranger-arc-mm10-2020-A-2.0.0, obtained from https://cf.10xgenomics.com/supp/cell-arc/refdata-cellranger-arc-mm10-2020-A-2.0.0.tar.gz). To capture CMO reads, the gene expression modality was separately processed using the cellranger count pipeline (v7.0.0 and v7.1.0), aligning it against the mm10 reference genome (refdata-gex-mm10-2020-A, obtained from https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz). Each nucleus was assigned to the sample of origin as part of further demultiplexing, for which a Gaussian mixture model was fitted to the CMO count distributions. The probability that a given CMO read count was derived from background or positive signal was then inferred (as previously described64). The resulting multiome data was further analysed in R (v4.3.0) or Python (v3.10.5) using Seurat65 (v4.3.0.9012) and Signac66 (v1.10.0). For reproducibility, the analysis was implemented as targets67 (v1.2.2) pipeline using renv (v1.0.2).
Pre-processing was done per 10X lane (that is, per multiplexing group of 4–6 samples). Peaks for the snATAC-seq modality were called for each run separately using MACS2 (v2.2.7.1) within Signac. A set of consensus peaks for the snATAC-seq analysis was created according to Signac. CMO demultiplexing information was taken from Cellranger output and added as metadata.
Doublets were identified using scDblFinder68 (v1.14.0), separately for snRNA-seq and snATAC-seq modality. snRNA-seq parameters for the scDblFinder function were dims = 50, cluster = TRUE and knownUse = ‘discard’. snATAC-seq features were first aggregated (aggregateFeatures function, k = 500) and parameters for scDblFinder were dims = 50 and aggegateFeatures = TRUE. Cells identified as doublets based on CMO data were included as known doublets (knownDoublets parameter). Quality control metrics such as the percentage of mitochondrial reads, nucleosome signal and transcription start site enrichment were calculated per run. After removing known and estimated doublets, outlier thresholds for snRNA-seq read counts (nCount_RNA) or snATAC-seq read counts (nCount_ATAC) were derived using the isOutlier function from the scuttle package69 (v1.10.2) with parameters nmads = 3, type = ‘both’ and log = TRUE. Nuclei falling outside these thresholds or with more than 5% of mitochondrial reads were removed.
Integration into combined analysis
Complete datasets for each time point or tissue were created by integrating individual Seurat objects representing separate 10X lanes (that is, multiplexing groups). snRNA-seq data were normalized using SCTransform70,71 from Seurat with parameter vst.flavor = ‘v2’. Integration of normalized snRNA-seq data was performed according to Seurat using the reciprocal principal component analysis (PCA) method with the following non-default parameters: SelectIntegrationFeatures(nfeatures = 2,500), FindIntegrationAnchors(reduction = ‘rpca’) and IntegrateData(dims = 1:50). Integration of snATAC-seq data was performed at the level of the LSI embedding using the IntegrateEmbeddings(dims.to.integrate = 1:50) function with the integration anchors derived from snRNA-seq data. Identification of cell types was performed using the combination of both modalities: snRNA-seq and snATAC-seq. A PCA embedding was derived for the integrated snRNA-seq data, the weighted nearest neighbour graph was calculated using the FindMultimodalNeighbors function with the integrated LSI and integrated PCA reductions as input and dims_atac = 2:50 and dims_rna = 1:50. The resulting weighted nearest neighbour graph was used for UMAP visualization of the dataset and identification of clusters using the FindClusters function with parameters algorithm = 3 (smart local moving algorithm), resolution = 0.5 (forebrain and cortex) or 0.3 (cerebellum). Cells designated as outliers per cell cluster were designated as such and removed using the Mahalanobis distance (more than 100) calculated from genes expressed in more than 80% of the cells.
Label transfer and cell-type annotation
To support cluster annotation, reference datasets for each tissue and time point were used (a list of respective references used is provided in Supplementary Table 1). Seurat objects were created from all datasets except where data were provided in this form. Normalization and dimensionality reduction were performed in alignment with the published studies where possible. Cell-type annotation of clusters was done in multiple steps. Label transfer from multiple published datasets (see above) was used for automated annotation of clusters using Seurat. Seurat’s FindMarkers was used to determine marker genes for each cluster. Cluster annotations were manually validated afterwards. Here information from different reference datasets, in addition to the expression of canonical marker genes, was validated. Clusters annotated as striatal cell types were excluded from further analyses due to sample preparation variabilities, that is, brain region adjacent to the forebrain, and the presence or absence thereof is highly variable. To reduce the complexity of downstream analyses, clusters belonging to the same cell-type annotation were merged. To avoid merging disparate clusters, differential expression (using FindMarkers of subgroups of clusters) was used as a reference.
Cell variability
Single-cell gene expression heterogeneity within each cell-type cluster and mutant was tested and visualized using coefficient of variation histogram plots and permutation testing. Permutation testing refers to the number of genes designated as significantly different between (1) cells of the same transgenic mouse line or the C57BL/6J WT mice, or (2) cells of different transgenic mouse lines or the C57BL/6J WT mice. Twelve cells were randomly sampled from the respective population (or populations) two times, and differences in expression values were obtained. Any difference value that was obtained from less than 10 cells with an expression value was discarded. Then, two such differences in expression value sets (either for (1) or (2) were compared with a t-test and deemed to have a high variability when their P value was below 0.05. This was repeated ten times. Furthermore, differences were also calculated between cell-type clusters of one transgenic mouse line and the C57BL/6J WT mice.
Differential abundance analysis
Differential compositional analysis of cell abundances between all mutants linked with ASD collectively and WT mice was conducted using the scCODA Python package72 (v0.1.8), using Hamiltonian Monte Carlo sampling. The reference cell types used were as follows: layer I (Cajal–Retzius) neurons for the cerebral cortex at E14.5, CGE interneurons (interneurons CGE (Vip)) for the cerebral cortex at P4, interneurons Sncg (Vip) for the cerebral cortex at P14 and choroid plexus cells for the cerebellum at P4. The credibility of the effects was evaluated at two FDR thresholds: 0.05 and 0.01. Differential abundance of cell-type clusters across individual mutants linked with ASD was also tested with edgeR73 (v3.42.4) using absolute cluster abundances in cells per sample as input. For cortical time points, the model included sex as a covariate.
Pseudotime analysis
To reconstruct the differentiation trajectory from RG to the cortical plate cells at E14.5, we created an embedding based on the transcriptomic profiles using the Python package Harmony74 (harmonypy, v0.0.5) for batch correction. After that, Palantir75 (v1.0.1) was used to generate a latent representation that preserves the pseudotime ordering of cells. Utilizing the first five components of this representation, the SimplePPT algorithm76,77 was applied to learn the principal tree that captures the differentiation process. Subsequently, the scFates tool78 (v1.0.6) was used to map this principal tree back onto the constructed UMAP using Harmony-corrected PCA coordinates. Pseudotime values for each cell were assigned by projecting the cells onto the principal tree. The RG cell cluster was defined as the root of the trajectory, progressing towards cells of the cortical plate. The pseudotime analysis included RG cells, subventricular zone migrating cells, intermediate progenitors and cortical plate cells.
Pseudobulk differential gene expression analysis
Differential gene expression between genotypes within cell-type clusters was performed on pseudobulk data. Either (1) all female mice only, (2) all male mice only, or (3) all female and male mice without sex indicators were used for testing between one transgenic mouse line and the C57BL/6J WT mice.
Pseudobulk profiles per cluster and per sample were derived using the Libra package (v1.0.0). Mitochondrial genes were removed, and only genes with at least 15 counts in at least two samples (single sex) or 30 counts in at least two samples (pooled sexes) were included. DESeq2 (ref. 79) (v1.40.2) was used to test for differential gene expression using test = ‘Wald’, fitType = ‘local’ and an adjusted P value threshold of 0.05. The model included sex as a covariate for the analysis of the forebrain and cortex. Shrunken fold changes for downstream rank-based analyses were calculated using the apeglm method80 implemented in the lfcShrink function from DESeq2. The tests involved three comparison types: (1) comparing each single mutant genotype to WT controls, (2) comparing each of the single mutant genotypes to all other mutant genotypes (excluding WT controls), and (3) comparing the WT mice to all mutants collectively. Composite differential expression analysis was carried out by comparing gene names designated as differently expressed (DESeq2) and organizing these as Venn diagrams. Furthermore, a leave-one-out cross-validation strategy, which systematically removed one of the C57BL/6J WT mice from the DESeq2 comparison, was used to investigate whether any of the WT mice were a strong outlier. To assess differences between female and male WT controls, WT controls were also compared with a similar statistical setup as described above using pseudobulk differential gene expression analysis. Here, instead of adjusted P values, nominal P values were used to assess the extent of the sex differences with a critical α-threshold of 0.05.
To evaluate the statistical significance of overlap in differential gene expression patterns between mutants, each mutant pair was compared against distinct sets of C57BL/6J control samples. Control samples (including both male and female samples) were randomly divided into two disjoint groups. For each mutant pair, differential expression analysis was performed in DESeq2 using the same model as described above (‘~group + sex’), except for Ptchd1 mutants, where all samples were male and the model ‘~group’ was used. Comparison of each mutant to a distinct control set ensured statistical independence between the two comparisons. Overlap between each mutant pair was assessed by comparing the top 300 most DEGs (ranked by P value), and the observed-to-expected ratio and statistical significance of the overlap were determined using Fisher’s exact test. Only genes with a minimal expression level (baseMean ≥ 1) in at least one mutant of a given pair were included in the analysis.
To compare the extent of transcriptional differences across mutants and cell types, we used a downsampling-based differential expression analysis aimed at balancing statistical power. Pseudobulk profiles were generated for each genotype–cell type combination and downsampled to 10,000 molecules per profile. Profiles with fewer than 95% of this target were discarded, and only combinations with at least five valid samples were included in the analysis. Differential expression testing was performed with DESeq2 using the design formula (‘~genotype + sex’). To mitigate stochastic effects from subsampling, we conducted 100 independent downsampling iterations, and for each mutant–cell type comparison, we reported the average number of genes with an adjusted P < 0.05.
Gene Ontology analysis and GSEA
Gene Ontology term analysis was performed using the enrichGO function from the clusterProfiler package81 (v4.8.2) using the list of DEGs using a cut-off P < 0.01 and cut-off q < 0.05.
GSEA was performed as previously described82, using the gseGO function from the clusterProfiler package81 (v4.8.2) with default parameters except minGSSize = 100 and ont = ‘ALL’ using a cut-off P < 0.05. All genes were used as input for clusterProfiler and the list of genes was sorted by log2FC using shrunk fold changes using the lfcShrink function. Enrichment was assessed using a Kolmogorov–Smirnov-like statistic with permutation-based significance testing. Benjamini–Hochberg FDR correction of P values was used for both analyses.
To remove redundancy in enriched gene ontologies, the simplify-method function from clusterProfiler was used to visualize enriched gene sets (Figs. 2b and 5b).
We performed Gene Ontology term enrichment analysis on genes associated with differential accessibility at the gene-body level. DAR (P < 0.01) per cell type were enriched with Gene Ontology terms using the clusterProfiler R package (using the Biological Process ontology, using a cut-off P < 0.01 and a cut-off q < 0.05. Enriched terms were illustrated using the enrichplot R package (v1.20.1) with the dotplot function.
Constructing enrichment maps using Cytoscape
To visualize the results of enriched gene sets and pathways for mutants linked with ASD, the EnrichmentMap plugin83 (v3.4.0) for Cytoscape (v3.10.1) was used. Enriched gene sets for excitatory neurons were loaded into the EnrichmentMap plugin with an FDR cut-off value of 0.1. Gene sets were not filtered, setting no minimum experimental cut-off (Extended Data Fig. 4).
Summarizing of Gene Ontology terms using REVIGO
REVIGO84 was used to consolidate GSEA terms into a reduced set of parent terms. For this, a comprehensive list of all GSEA terms, irrespective of time point or genotype, was constructed and uploaded. The parameters ‘small (0.5)’ and ‘remove obsolete GO terms’ were used, resulting in robust and comprehensive clustering. Subsequently, all detected GSEA terms were mapped to their new parents. The parent terms were organized in subclusters of circular heatmaps, which were implemented with the circlize (v0.4.15) R package85. Representative terms were highlighted in the corresponding figures showing circular heatmaps.
Synaptic function and localization analysis using SynGO
Synaptic Gene Ontology (SynGO)86 (v1.2) was used to identify types of known synaptic proteins (in presynaptic and postsynaptic compartments) and different functional classes. Lists of DEGs of the respective mutants linked with ASD of P4 and P14 excitatory neurons were used as input to investigate overrepresented synaptic terms. Biological processes were selected for enrichment.
Differential accessibility
We conducted pseudobulk testing using the DESeq2 R package79 (v1.34.0) to identify peaks that differ between ASD mutant and WT mice. Given the low signal-to-noise ratio in scATAC-seq experiments87 and the high number of comparisons, we used an uncorrected P value threshold of 0.001 to include peaks in the list of DARs. Motif enrichment analysis was performed on a set of downregulated peaks using the SnapATAC2 Python package (v2.3.0), with binomial testing and Benjamini–Hochberg FDR correction of P values. We tested for correlations between the gene-body differential accessibility P values and the corresponding gene lengths. Our analysis revealed no significant association (<0.1), demonstrating that gene length is not a systematic confounding factor.
Network-based analyses
Generation of reference protein–protein interaction network
Several sources were merged to generate the protein–protein interaction network (PPI) used in this study. In brief, interaction data for mouse was downloaded from MIPPIE88, IntAct89, STRING90, Reactome91 and BioGRID92. For STRING, only physical interactions were used. The BioGRID database was downloaded separately, as newer releases were not integrated in the MIPPIE database. In addition, we incorporated mouse-brain-specific interaction data from Skinnider et al.93 to ensure coverage of target tissue-specific interactions. The resulting PPI comprises 18,921 proteins with 417,767 undirected interactions.
Network colocalization of DEGs
Raw count data were log normalized using a library size normalization \(\log (\frac{x}{\mathrm{libary}}\times 10,000+1)\). DEGs from pseudobulk differential gene expression analysis of all 11 mutant lines were evaluated for similarity in their underlying PPIs using a modified version of a previously published method94 applied to our reference PPI network. Key modifications include the use of a symmetric random walk, which inherently corrects for degree distribution95, and an extension of the original tied diffusion algorithm96 to handle more than two random walks. In brief, for each genotype and time point, we performed symmetric random walks using the pseudobulk DEGs as seed nodes. For each walk, network nodes were ranked in descending order according to their visiting probability. An aggregated rank for each node was then computed as the geometric mean of its ranks across all walks, referred to as the rank product97. Only nodes with a rank product below 3,000 were retained for further analysis. This approach builds on NetColoc, originally designed to compare pairs of diseases using the tied diffusion algorithm introduced by Paull et al.96 to identify shared components of molecular networks between two sets of molecules. Although NetColoc was only demonstrated for pairwise comparisons, we argue that the method is not inherently restricted to two sets, as described below. In brief, TieDIE begins by performing a random walk on a common reference interaction network for each set of molecules of interest. Each random walk assigns weights to network nodes according to their proximity to the respective molecules. The shared neighbourhood among the resulting walks is then determined by applying a threshold to a weighting function \(z=f(R)\), where R denotes the set of random walk outputs and \(f(.)\) is an arbitrary function that prioritizes nodes with consistently high weights across walks. To accommodate an arbitrary number of random walks, we adopted a previously published ranking approach developed to identify highly ranked DEGs across multiple differential expression analyses97. Specifically, our ranking function is the rank product across all random walks in R, calculated as the geometric mean of the rank for each node across the walks.
Computing gene systems from colocalized genes
Colocalized gene sets were grouped into clusters of genes with similar molecular functions, referred to as gene systems, using HiDeF98. We applied HiDeF to cluster colocalized genes separately for each time point. To focus on more specific systems, we excluded clusters containing more than 250 genes, as such large clusters are often nonspecific artefacts of the hierarchical clustering used by HiDeF. We further removed clusters with a persistence score below 20 or without any overlap with the corresponding seed gene sets. For the remaining clusters, we assessed functional enrichment in GO_BP_2023 (Biological Processes), GO_MF_2023 (Molecular Function), GO_CC_2023 (Cellular Component) and MGI_Mammalian_Phenotype_Level_4_2024 using the ENRICHR API via gseapy99. To avoid redundant gene systems across time points, we merged clusters based on their enrichment profiles using hierarchical clustering. The similarity between two clusters was calculated as the jaccard index of the tokenized enrichment terms. If no enrichment was available for either cluster, similarity was computed from the overlap with genes in the system.
Gene system–genotype associations
Associations between gene systems and genotypes were assessed using a regression-based method described in Jin et al.100. In brief, we first filtered the gene systems to retain only those with at least four genes present in our snRNA-seq dataset, removing any system that fell below this threshold after excluding genes absent from the data. For each time point, we then assessed the association between gene systems and genotypes by fitting regression models of the genotype variables on the module eigengene scores of each system, calculated per cell in the corresponding subset of the snRNA-seq data. The full model was fit in statsmodels (v0.14.3) using the WT as reference, and is specified as follows:
module eigengene score ~ genotype + nFeature_RNA
The resulting regression coefficients represent the associations between each genotype and gene system relative to WT. To further reduce the number of systems for visualization, we retained only those that (1) were associated with more than 20% of genotype–time point combinations, and (2) exhibited substantial variation across genotypes and time points, defined as having a maximum absolute difference from the median greater than 0.2 across all combinations (genotype + time point).
Network overlap computation
We computed pairwise gene set overlaps on the PPI as introduced in Menche et al.101. Gene sets were defined for each genotype and time point by ranking genes within the PPI network according to their absolute log2FC values from pseudobulk differential expression analysis. Excitatory and inhibitory neurons were analysed separately. Network overlap analyses were then performed for all pairs of genotype–time point combinations following Menche et al.101. For each comparison, overlaps of 1,000 randomly drawn gene set pairs were computed to generate a reference distribution. To ensure comparability, we implemented a degree-matched randomization procedure by binning PPI nodes into 10 equally sized bins according to their degree and drawing random genes from these bins based on the degree distribution of the actual gene sets. Statistical significance of observed overlaps was assessed by computing empirical P values, defined as the proportion of random overlaps more extreme than the observed overlap. Resulting P values were corrected for multiple testing using the Benjamini–Hochberg procedure as implemented in the statsmodels package102 (v0.14.3). Coefficients different from zero were assessed using two-sided t-tests with Benjamini–Hochberg correction. Coefficients with adjusted P < 0.1 were considered significant and indicated with a black outline in the plot.
Bulk RNA-seq
Sample preparation
Brain tissue from Bckdk−/−, Cul3+/−, HnrnpU+/− and C57BL/6J colony-matched controls (n = 3 per sex + genotype) across time points (E14.5, P4 and P14), as well as from E16.5 HnrnpU and WT littermate mice (n = 3 per sex + genotype) were used for the bulk RNA-seq experiment. Animals coming from at least three separate litters were collected to avoid potential litter-specific biases. Timed-pregnant females were decapitated, and embryos were rapidly dissected on ice. Cortical tissue was dissected in ice-cold PBS, for embryonic tissue, the meninges were removed, and the tissue (whole brain for E16.5 HnrnpU animals, forebrain for mutants linked with ASD at E14.5, P4 and P14) was snap-frozen in liquid nitrogen and stored at −80 °C. RNA extraction for the E16.5 HnrnpU dataset was performed using TRIzol–chloroform extractions as described in the Methods section for RNA isolation. For the dataset using multiple mutants linked with ASD at E14.5, P4 and P14, the RNAeasy Mini RNA isolation kit (74104, Qiagen) was used according to the manufacturer’s protocol. One microgram of RNA obtained was used for cDNA library preparation using the NEBNext Ultra II Directional RNA Library Prep kit (E7765, NEB) with unique dual indexing primers (NEBNext Multiplex Oligos for Illumina, E6440S, NEB) and poly(A) selection with the NEBNext Poly(A) mRNA Magnetic Isolation Module (E7490, NEB), according to the manufacturer’s protocol. RNA quality, library quality and size were assessed with the RNA 6000 Nano Kit (5067-1511, Agilent) and High Sensitivity DNA Analysis Kit (5067-4626, Agilent) and the Bioanalyzer 2100. Libraries were sequenced on the Illumina NextSeq 2000 and NovaSeq X platforms by the Vienna Bio Center NGS facility.
Pre-processing and analysis of RNA-seq
For the analysis of E16.5 HnrnpU mutant and WT animals, demultiplexed FASTQ files were filtered and trimmed using Cutadapt103 (v4.1) and aligned to the mm10 reference genome using STAR104. Similarly, this method was used to obtain read counts. All further downstream analyses were performed in R. The same pipeline used for pseudobulk DGE analysis of the snRNA dataset as described above was used. However, the pseudobulking step was omitted, and only the mutant genotype was compared with WT littermate mice. One sample (HET genotype, female) was removed due to an increased number (approximately 10%) of duplicated sequences. For the analysis of mutants linked with ASD (Bckdk−/−, Cul3+/− and HnrnpU+/−) and C57BL/6J colony-matched controls, data were analysed using R (v4.5.0). Demultiplexed FASTQ files were filtered and trimmed using Cutadapt103 (v4.1). Samples were aligned to the GRCm39 reference genome using STAR (v2.7.11b), and gene-level counts were obtained using FeatureCounts (v2.1.1)105. Mutant samples were compared with WT C57BL/6J colony-matched controls. Unwanted variation was corrected using the RUVSeq package (v1.42.0)106. To define empirical control genes, an initial DESeq2 analysis was performed, and all genes were ranked by P value. The top 5,000 DEGs were excluded, and the remaining genes were used as empirical controls for RUV (remove unwanted variation) correction. The value of k was selected as the lowest number that effectively reduced unwanted variation, based on visual inspection of correction effects in the PCA space, following the recommendations of Risso et al.106.
RNA isolation and quantitative real-time PCR
Brain tissue from male and female Ash1l+/−, HnrnpU+/−, Kdm6b+/−, Kmt5b+/−, Trip12+/−, Usp7+/− and Wac+/− and their WT littermates were used to validate the expression of targeted genes. For RNA extraction of cortical samples, 700 µl TRIzol (Invitrogen) and 150 µl chloroform (Sigma) were used for homogenization, followed by centrifugation at 12,000g for 15 min at 4 °C. For cerebellar samples, 600 µl TRIzol and 130 µl chloroform were used.
The upper aqueous phase was transferred to a new tube, and 1.5 volumes of 100% ethanol (EtOH) were added. RNA was extracted using Zymo-Spin IC columns (Zymo Research). The aqueous phase–ethanol mixture was loaded onto the column and washed with 400 µl of 70% EtOH following treatment with RQ1 DNAseI (Promega; 5 µl + 5 µl reaction buffer + 40 µl 70% EtOH) for 15 min at room temperature. After two washes with 70% EtOH, the sample was eluted from the column with DEPC-treated H2O. RNA concentration was measured using the NanoDrop spectrophotometer (Thermo Scientific). Of RNA, 1 µg was used for cDNA preparation with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher). cDNA was diluted 1:5 and used for quantitative PCR analysis with Lightcycler 480 Sybr green master mix (Roche) on a Real-Time PCR Roche Lightcycler 480 machine. Samples were run in triplicates for target genes, and a list of all oligonucleotides used is provided in Supplementary Table 1. ΔCq expression levels (relative mRNA) were calculated upon normalization to housekeeper genes and plotted.
Cytometry, imaging and analysis
Sample preparation and imaging of EdU-injected and CUBIC-cleared brains
Timed-pregnant female mice were injected with 0.05 mg g−1 of 5-ethynyl-2′-deoxyuridine (EdU) at E14.5 and euthanized 6 h later; embryos were decapitated, rapidly dissected on ice and immersion fixed in 4% paraformaldehyde overnight at 4 °C, followed by a washing step in 1X PBS. Cells in S-phase that incorporated EdU were detected using click chemistry (BCK-EdU555IM100+IV-M, Baseclick). Embryonic heads were cleared using the CUBIC clearing protocol107 with adaptations coupled with EdU staining. In brief, fixed, washed samples were incubated for 4 h at 37 °C on a rotating shaker in CUBIC-L/R1a in mH2O. Samples were incubated in CUBIC-L/R1a solution for 48 h at 37 °C on a rotating shaker with a change of solution after 24 h. Afterwards, samples were washed with PBS and incubated overnight at 37 °C with PBS–0.2% Triton X-100–20% DMSO–0.3 M glycine. Before incubating the samples with the EdU reaction solution according to the manufacturer’s protocol overnight, samples were incubated for 2 h in PBS–0.2% Triton X-100 and briefly washed twice with PBS. After incubating with the EdU reaction, samples were again washed twice with PBS and 2 h in PBS–0.2% Triton X-100. To match the refractive index, samples were stepwise incubated for 2 h at 37 °C in ascending concentrations of CUBIC-R + (N) in mH2O. In the case that embryonic heads were not fully transparent, they were incubated overnight in CUBIC-L/R1a at 37 °C on a rotating shaker, followed by stepwise refractive index (RI) matching. Cleared heads were imaged in oil (Gelest #PDM-7040, RI of 1.556; and J62592.AP, Thermo Fisher, RI of 1.467). Both oils were mixed to match the refractive index (RI ~ 1.52) of the cleared heads. Imaging was performed on a Zeiss Lightsheet 7 microscope equipped with dual-side illumination and two camera detection (PCO edge 4.2 sCMOS) using ZEN Black imaging software (v3.10) using 5×/0.1foc illumination and 5×/0.16 EC-Plan Neofluar detection objectives with the following settings: filter set: band pass 505-530/LP560/LP585, Zoom 1.4, lightsheet thickness of 7.76 µm, scaling xyz of 0.67 µm × 0.67 µm × 1.5 µm. Imaging was performed using single-side illumination to achieve optimal sample illumination and to prevent potential duplication of nuclei in the images. Whenever possible, both the left and right brain hemispheres were imaged, and damaged or incompletely cleared samples were excluded from the analysis. Images of EdU-stained, cleared embryo heads were analysed automatically using Zeiss Arivis Pro software (v4.2.0). A semi-automated analysis pipeline was optimized to count EdU-positive cells within one cortical hemisphere, defined as a region of interest, using the Blob Finder function implemented in Zeiss Arivis Pro (v4.2.0).
Sample preparation and analysis of EdU-injected brains by cytometry
Timed-pregnant female mice were injected with 0.05 mg g−1 EdU at E14.5 and euthanized 2 h later; embryos were rapidly dissected on ice. Cortical tissue was isolated in ice-cold HBSS (14175053, Thermo Fisher), meninges were removed and forebrain structures were dissociated using an established papain-DNAse protocol108. Following dissociation, S-phase cells that incorporated EdU were detected by click chemistry (BCK-EdU488FC50+IV-S, Baseclick). In brief, washed samples were fixed for 15 min at room temperature and permeabilized using the supplied 1X saponin-based permeabilization and wash buffer. Two washes (1% BSA in PBS, 300g for 10 min each) were performed between fixation and permeabilization. Following fixation and permeabilization, the click reaction was prepared according to the manufacturer’s protocol, and cells were incubated for 30 min at room temperature, protected from light. Samples were washed twice with supplied 1X saponin-based permeabilization and wash buffer and centrifuged at 300g for 10 min per wash. A control sample was processed identically up to the permeabilization step and used as negative control for cytometry analysis. Samples were analysed on a Beckman Coulter Cytoflex LX. Cell populations were gated based on forward scatter versus side scatter, and doublets were excluded via forward scatter area versus height gating. Fluorescence of EdU in gated single cells was measured in the B525-FITC-A channel (excitation of 488 nm and emission of 525/40 band pass. Total events (n = 50,000) were acquired per sample, and gating and analysis were performed using Beckman Coulter’s CytExpert software (v2.7).
RNAscope assay
Spatial gene expression analysis was performed using the RNAscope94 Multiplex Fluorescent v2 Assay (323110, ACD) kit including specific probes targeting Pdgfra and Slc1a3 mRNA (480661 and 430781-C3, ACD). For sample preparation, P4 mice were decapitated, and the brain was dissected rapidly on ice. The cerebellum was removed, and the brain was embedded in pre-cooled Tissue-Tek OCT Compound (Sakura Finetek) and stored at −80 °C until further use. Tissue was sliced at 20 μm at a CryoStar NX70 cryostat (Thermo Scientific) and directly mounted on SuperFrost Plus Adhesion slides (Epredia). Sections were stored at −80 °C until use. The assay was performed according to the instructions provided by the RNAscope Multiplex Fluorescent v2 Assay kit. In brief, sections were fixed in 4 °C cold 4% paraformaldehyde and further pre-treated using a H2O2 and protease digestion. After pre-treatment, sections were incubated with the probes targeting the mRNAs of interest. Probes were further hybridized to a cascade of signal amplification molecules, followed by a hybridization with TSA vivid fluorophore 520 for Pdgfra and TSA vivid fluorophore 570 for Slc1a3 (323271 and 323272, ACD). After the hybridization steps, sections were stained with a nuclear counterstain and mounted using ProLong Gold Antifade Mountant (P36930, Invitrogen). Images were acquired using a LSM900 inverted confocal microscope using 10× objective and analysed in ImageJ. To assess the number of astrocytes and OPCs, Slc1a3+ and Pdgfra+ cells were manually quantified in blinded images of coronal brain sections of similar size in at least three images per mouse.
Electrophysiology
Brain slices were obtained from P7 Bckdk−/−, HnrnpU+/−, Trip12+/−, Usp7+/−, littermate and colony-matched C57BL/6J control mice. Acute coronal slices (300 µm) were prepared from the primary somatosensory cortex. Animals were decapitated and whole brains were rapidly removed from the skull and sectioned using a VT 1200S vibratome (Leica Microsystems). Slices were then allowed to recover in in regular artificial cerebrospinal fluid in 1 mM Ca2+ at 35 °C for 10 min, followed by a cool-down to room temperature for at least 20 min in artificial cerebrospinal fluid, containing: 118 mM NaCl, 2.5 mM KCl, 1.5 mM MgSO4 × 7H2O, 1.25 mM NaH2PO4, 2 mM CaCl2, 20 mM sucrose, 3 mM myo-inositol, 10 mM glucose and 26 mM NaHCO3 (320 mOsm l−1, pH 7.35–7.40).
Patch pipettes (3–5 MΩ; Sutter Instrument) were pulled on a P-1000 puller (Sutter Instrument) and filled with the intracellular recording solution, containing: 120 mM K-gluconate, 20 mM KCl, 2 mM MgCl2, 10 mM HEPES, 0.5 mM EGTA, 2 mM MgATP, 0.2 mM NaGTP, 23.4 mM sucrose and 0.3% w/v biocyin; adjusted to a pH of 7.35, and osmolarity adjusted to 302 mOsm l−1.
Excitability (current versus action potential frequency) was measured in current clamp mode by applying current steps from −10 to +50 pA in 5-pA increments (1-s duration and 8-s inter-sweep interval). Input resistance (Rin) was calculated as the linear regression of the current–voltage relationship from −10 to 0 pA. Membrane time constant tau was measured by an exponential fit of the initial voltage change in response to a −10-pA current injection. Membrane capacitance was calculated by dividing tau by Rin. Resting membrane potential was measured as the average voltage at the 0-pA current injection sweep. Action potential properties were measured in the first evoked action potential, with threshold potential measured at the onset of the steep rising phase. Action potential amplitude was measured as the difference between threshold and peak potential. Action potential full-width at half-maximum was measured as the time between rising and falling phase at half of the maximal amplitude. Action potential afterhyperpolarization was measured as the voltage 50 ms after the peak relative to the threshold potential. mEPSCs and mIPSCs were recorded for 3–5 min in voltage clamp mode at a holding potential of −60 mV and in the presence of 1 µM tetrodotoxin + 10 µM bicuculline methiodide (mEPSCs) or 1 µM tetrodotoxin + 10 µM DNQX + 50 µM D-AP5 (for mIPSCs). Access resistance (Ra) was monitored throughout all recordings, and recordings with Ra exceeding 20 MΩ were excluded. Miniature neurotransmission was manually analysed using Mini Analysis Program (Synaptosoft, v6.0.7). Passive and active membrane properties were analysed using Clampfit (v11.3.0.02). All recordings and analyses were performed with the experimenter blind to the animal experimental condition.
Minimal conductance-based model
To model the electrophysiological properties of WT neurons, we utilized the biophysical model of a layer II/III pyramidal neuron from the mouse visual cortex provided by the Allen Brain Institute (https://modeldb.science/184162). Membrane properties and morphology were modified to match the developmental stage of P7 and electrophysiological characteristics (Rin and membrane time constant) observed in our data. The soma and dendrites were scaled to 70% of their original radius and length, and the axon was reduced to 50% of its original dimensions. To match the observed frequency–current relationship, the membrane capacitance was set to 0.38× and the passive membrane conductance to 0.5× of the original values in their respective compartments. Reversal potentials were set to −75 mV for passive current, +30 mV for Na+ and −110 mV for K+ ions. Simulations were performed in Python 3 using the NEURON simulator109. Following the described parameter adjustments, the model exhibited a net input resistance of approximately 988 MΩ, a membrane time constant of 63 ms and a resting membrane potential of −61 mV, approximately matching the mean passive properties of a WT control neuron.
The model includes multiple ion channel types: persistent potassium (Kp), delayed rectifier potassium (Kt), A-type potassium (Kv3.1), Ca2+-activated potassium (SK), non-inactivating potassium (M-channel), persistent sodium, transient sodium, hyperpolarization-activated cyclic nucleotide-gated (Ih), Im channels, high-voltage-activated calcium and low-voltage-activated calcium channels.
Behavioural analysis of Bckdk-mutant and Kmt5b-mutant animals
Gait recordings
Bckdk-mutant and WT littermate animals (P13–15) were acclimated to the testing room 1 h before behavioural recording and were not handled before testing due to their young age. Animals were led to explore an empty open-field arena (50 × 50 × 40 cm) constructed from IR-transparent acrylic covered with a lid, creating a dark environment during testing. Animals were recorded from below using an IR-camera (Basler acA3088-57µm) ensuring continuous visibility of the paw of the animals. Video acquisition was performed with OBS Studio (v31.0.1). Each recording lasted a minimum of 15 min to capture a sufficient number of stride cycles. The experimenter was blind to the genotype during recordings. Pairs of control-mutant littermates were randomly tested in the morning and in the afternoon to control for circadian rhythm.
Analysis of mouse pose estimation and gait
Body parts of mice were tracked using SLEAP110 (v1.4.1) by defining a SLEAP skeleton containing 14 body parts. Four are at the head of the mouse (nose, ear left, ear right and neck base), the four paws (forepaw left, forepaw right, hindpaw left and hindpaw right), at the centre and sides of the mouse body (spine centre, left side, right side) and three parts for the tail (tail base, tail middle and tail tip).
After annotating 364 frames of a single mouse across WT and mutant conditions, we trained both a centroid and a centred instance model using the top-down pipeline in SLEAP. In the centroid model, which detects the central body region, we used an input scaling of 0.5, an output stride of 2, and defined spine centre as the anchor point with σ = 2.5 px. All other network parameters were kept at their default values. The model was trained for 128 epochs with early stopping enabled (plateau patience = 32 epochs).
The centred instance model, which detects all body parts within a bounding box centred on the anchor point predicted by the centroid model, was configured with an input scaling of 1.0, output stride of 2, maximum stride of 32, number of filters of 32 and the filter rate of 2. We used a σ = 1.75 for the body parts. This network was trained for 500 epochs allowing for early stopping (plateau patience = 84 epochs).
Both models used a UNet backbone and were trained with an Adam optimizer (initial learning rate = 0.0001). Data augmentation included random rotations, horizontal flips and random size variations between a scale of 0.95 and 1.05. The validation dataset comprised 10% of the total annotated frames. Training was performed on a NVIDIA RTX 4080 SUPER GPU with 16 GB of VRAM.
The centroid model achieved a detection precision and recall of 1 on the validation set. The centred instance model reached a body-part recall of 0.992 and precision of 0.996, with an average pixel error of 2.59 px. To track the detected mouse poses, the default centroid tracker in SLEAP using the Hungarian method for instance matching in a frame window of size 5 was used.
To analyse mouse gait, we implemented mouse gait and posture metrics as previously described111. To extract strides, we first egocentrically aligned each mouse pose by centring the mouse to its central part (that is, spine centre) and rotating the coordinate system such that the neck base pointed upwards along the y axis. In this egocentrically aligned coordinate system, local maxima of the y location of paws (for example, forepaw left) were identified as the time point of foot-strike events, marking the start of the stance phase. Similarly, local minima determine the toe-off event, marking the beginning of the swing phase.
Before detecting extrema, paw trajectories were smoothed with a one-dimensional Gaussian filter (σ = 3 frames, corresponding to 0.05 s). Maxima and minima were then extracted using the scipy.signal.find_peaks function, requiring a minimum peak prominence of 11 px (0.6 cm).
Strides extracted from one paw part were validated by the corresponding opposite paw. A stride was considered valid if, during its stance phase, the opposite paw exhibited a foot-strike event. Both forepaw and hindpaw strides were validated and included in downstream stride analyses.
To focus on periods of active locomotion, valid strides were filtered to include only those occurring during forwards movement. Forwards motion was quantified by projecting the displacement vector (instantaneous speed) of the mouse onto its body axis (defined by the unit vector from spine centre to neck base). This resulting forwards speed was smoothed using a Gaussian kernel (σ = 60 frames; 1 s) and thresholded with 0.5 cm s−1 to define periods of forwards movement. Only strides occurring during forwards movement were used for gait metric computation. All raw stride, step and speed metrics were normalized to the body length of the mouse, defined by the distance from tail base to nose. All resulting stride metrics were averaged per movie.
Stride extraction and metric computation were implemented using the TadPose library112 (v0.4.0), a Python library for the analysis of SLEAP pose-estimation data. Mouse forwards movement was calculated in custom Python Jupyter notebook using the scipy library (v0.15.3). Left and right forepaws and hindpaws were analysed separately but averaged to obtain overall gait metrics distinguishing between forepaw and hindpaw movements.
Statistics
Information on the statistical analyses is provided in the corresponding Methods sections. Blinding and randomization were not applied unless otherwise specified. Each n value refers to the number of biological replicates (that is, animals or cells) used per condition.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
