Transcriptomic convergence and the female protective effect in autism

KD with CasRx

CasRx¹⁷, specifically NLS-RfxCas13d-NLS (Addgene, 109049), was used in combination with target-specific gRNAs for KD of targeted transcripts. A deep-learning-based tool⁶⁶ was used to design gRNA spacer sequences targeting the genes of interest. These spacer sequences were filtered to exclude homopolymers (more than four As, more than three Cs, Gs or Ts). Spacer sequences with potential off-target interactions were identified using BLAST⁶⁷ and excluded. For each target gene, three separate spacer sequences were selected and combined into a tandem pre-gRNA array, with each spacer sequence preceded by the CasRx direct repeat (DR36) sequence⁶⁶: CAAGTAAACCCCTACCAACTGGTCGGGGTTTGAAAC. Tandem pre-gRNA arrays were processed by CasRx into individual gRNAs for KD¹⁷. See Supplementary Table 1 for full sequences.

Molecular cloning

To create constructs for the initial NPC KD experiments, we used the lentiviral pSicoR plasmid⁶⁸ (Addgene, 11579) as the backbone. The mouse U6 promoter was replaced with the human U6 promoter, which was used to drive expression of the CasRx gRNA arrays. These constructs were used in combination with a CasRx knock-in cell line (described below).

To obtain enough NPCs for CasRx KD validation by auto-western blot (described below), an all-in-one construct was created that contained CasRx and EGFP expressed from the CAG promoter element and gRNA arrays expressed from the human U6 promoter. The backbone also contained inverted terminal repeats, which allowed the construct to be integrated into the genome by piggyBac transposase⁶⁹. These constructs were based on the LiOn design⁷⁰ (Addgene, 154019), such that expression of CasRx and EGFP was dependent on recombination of the construct by piggyBac transposase, preventing expression from the episomal plasmid. This all-in-one construct was also used for the fluorophore KD assays and NPC proliferation assays. For expression of fluorophores in NPCs, iOn-based constructs were created with mKO2, mScarlet or iRFP670 driven by the CAG promoter element.

For shRNA-based KD, shRNAs were designed using the Horizon siDesign Center tool (https://horizondiscovery.com/en/ordering-and-calculation-tools/sidesign-center). Predicted shRNAs starting with G were selected, and those with potential off-target interactions were identified using BLAST⁶⁷ and excluded. The selected shRNAs (Supplementary Table 13) were cloned into the pSicoR plasmid⁶⁸ (Addgene, 11579).

For organoid experiments, we created the FLICK construct by engineering several features into a backbone incorporating the FLEx design⁷¹. In particular, we included a modified Cre recombinase that included an intron, which prevented recombination of the plasmid during cloning⁷². We also included CasRx¹⁷ as well as all necessary components of the Tet-On 3G system⁷³ to enable doxycycline-inducible recombination. Furthermore, we included identifying barcode sequences a short distance upstream of the polyA sites for the fluorescent reporters, which enabled us to determine which KD construct was incorporated in each cell through scRNA-seq. Finally, we incorporated inverted terminal repeats, allowing the construct to be integrated into the genome by piggyBac transposase⁶⁹. The sequences of all constructs were verified by whole-plasmid sequencing (Plasmidsaurus).

Cell culture

All human cell culture work was performed under protocols approved by Boston Children’s Hospital IRB (S07-02-0087). The control male human iPS cell line 280 was previously generated and fully characterized⁷⁴. Human iPS cells and iPS cell-derived NPCs were grown on tissue culture plates precoated with Matrigel (Corning, 354277) or Geltrex (Thermo Fisher, A1413302). The cells were passaged as needed using Accutase (Millipore Sigma, SCR005) or ReleSR (STEMCELL Technologies, 100-0483). Human iPS cells were cultured in mTeSR Plus medium (STEMCELL Technologies, 100-0276). Human NPCs⁷⁵ were generated from iPS cells using a STEMdiff SMADi Neural Induction Kit (STEMCELL Technologies, 08582) by following the monolayer culture protocol for at least 3 weeks. NPCs were continuously grown in neural induction medium plus SMADi and were not switched to neural progenitor medium. The medium was supplemented with 10 μM Y-27632 (STEMCELL Technologies, 72304) for approximately 1 day following passaging of iPS cells or NPCs.

Cerebral organoids were generated from 280-iPS cells using a STEMdiff Cerebral Organoid Kit (STEMCELL Technologies, 08570), following some recommendations from Giandomenico et al.⁷⁶ with further modifications. In particular, the organoids were embedded in Matrigel (Corning, 354277) at day 7 and grown in expansion medium until day 14, at which point they were removed from Matrigel and switched to maturation medium on an orbital shaker. Full details are available at protocols.io (https://doi.org/10.17504/protocols.io.x54v95p5ql3e/v1).

None of the cell lines was authenticated. All cell lines routinely tested negative for mycoplasma contamination.

Derivation of CasRx knock-in iPS cells

CasRx¹⁷ was knocked into the AAVS1 safe-harbour locus of 280-iPS cells using targeting constructs that have been previously described⁷³. Individual iPS cell colonies were picked and expanded. Southern blot analysis was performed by Celplor to confirm correct targeting of the AAVS1 locus. Karyotype analysis was performed using KaryoStat (Thermo Fisher); this revealed a 2.4-Mb gain on chromosome 1, described in microarray nomenclature as arr[GRCh37] 1q32.1(204108163_206539638)x3. This region did not contain any of the target genes.

Lentivirus

Lentiviruses were generated using standard protocols that have previously been described⁷⁷. Lentiviruses were applied at a multiplicity of infection of 4.0, as this value yields more than 98% transduced cells. RNA was extracted 24 h later.

Library preparation for RNA-seq

Bulk RNA was extracted using a Direct-zol RNA Microprep Kit (Zymo Research, R2062), and 3′ RNA-seq libraries were prepared following a modified version of the TagSeq protocol⁷⁸, incorporating recommendations and pooling strategies from the CheapSeq protocol⁷⁹. Briefly, 200 ng of total RNA from each sample was combined with an indexed oligo(dT) primer and buffer (NEB, M0466L) as described for CheapSeq. RNA was fragmented in a preheated thermocycler at 94 °C for 2.5 min without the lid, after which the samples were moved immediately to ice. Each sample was then mixed with a template-switching oligo that contained an index and a unique molecular identifier (UMI), as well as other reagents required for reverse transcription, as described for CheapSeq. The reverse transcription reaction was performed at 42 °C for 90 min, without heat inactivation. Remaining single-stranded RNA and oligos were then digested by addition of 0.5 μl of buffer EB (Qiagen, 19086), 0.5 μl of thermolabile ExoI (NEB, M0568S) and 0.25 μl of RNase A (Thermo Fisher, EN0531) and incubation of the samples at 37 °C for 30 min, followed by heat inactivation at 85 °C for 5 min and a hold step at 4 °C.

Up to 8 samples with different index sequences were then pooled, and cDNA was purified by performing a 1.4× left-sided solid-phase reversible immobilization (SPRI) selection (Beckman Coulter, B23317). Q5 UltraII Mastermix (NEB, M0544L) was then used, with partial read 1 and partial read 2 primers, to amplify the cDNA for 7 cycles on a C1000 Touch Thermal Cycler (Bio-Rad, 1851197), following the CheapSeq PCR parameters. Amplified cDNA was purified through a 0.7× left-sided SPRI selection (Beckman Coulter, B23317). Final library construction was performed using Q5 UltraII Mastermix (NEB, M0544L) with standard TruSeq indexing primers (NEB, E6440S), with amplification for eight cycles following the CheapSeq PCR parameters. Libraries were purified by 0.4× to 0.65× double-sided SPRI selection (Beckman Coulter, B23317). Purified libraries were assessed using a TapeStation HSD1000 (Agilent, 5067-5584 and 5067-5585). For oligo sequences, see Supplementary Table 14.

The CheapSeq PCR parameters were as follows: 98 °C for 30 s; X cycles (see above) of 98 °C for 10 s, 65 °C for 75 s; 65 °C for 5 min; and 4 °C hold.

DEG analysis of bulk RNA-seq data from CasRx KD NPCs

Three distinct NPC samples were analysed per condition. Individual KD conditions were compared with three NTC conditions, which were considered together as a set of nine control samples. Gene expression was quantified using alevin-fry with UMI deduplication (salmon-v1.9). A salmon index and transcript-to-gene map were built using GENCODE v.45. Counts from distinct transcripts that mapped to the same gene were summed. Only genes with a baseMean expression greater than 10 in NTC samples were considered for further analyses. Differential expression analysis was performed using DESeq2 (v.1.42.0), with the model ~Batch + sgRNA.

For KD confirmation, P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing for the set of 20 directly targeted genes. The raw P value for SRCAP was derived from droplet digital PCR (ddPCR), as SRCAP was filtered from the RNA-seq analysis owing to low expression. See Supplementary Table 2 for details. For transcriptome-wide analyses, P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing for the set of 11,168 analysed genes (baseMean > 10 in NTCs).

Droplet digital PCR

To confirm KD of SRCAP, the RNA samples that had been used for RNA-seq (above) were assessed by ddPCR. RNA was first converted to cDNA using SuperScript IV VILO Master Mix (Thermo Fisher Scientific, 11756050). Next, we performed ddPCR by combining cDNA at an RNA-equivalent concentration of 0.25 ng μl⁻¹ (5 ng per 20-μl reaction), primers at a final concentration of 0.1 μM and QX200 ddPCR EvaGreen Supermix (Bio-Rad, 1864033). The ddPCR reactions were performed using a QX200 AutoDG Droplet Digital PCR System (Bio-Rad). Expression levels were normalized using the expression of housekeeping gene HPRT1. The results were analysed using the model Expr ~ Batch + type, where type was either NTC or the targeted gene (SRCAP). The resulting raw P value was adjusted with those for the other KD targets (from RNA-seq) using Benjamini–Hochberg correction for multiple hypothesis testing for the set of 20 directly targeted genes. See Supplementary Table 2 for details.

Primer design

Primers (Supplementary Table 15) were designed with Primer-BLAST⁸⁰, which assesses potential off-target interactions using BLAST⁶⁷. Oligos were ordered through Azenta or IDT. Primer specificity was confirmed by Sanger sequencing of PCR products (Azenta).

Auto-western blot

For analysis of protein expression levels following CasRx KD, LiOn-based KD constructs (‘Molecular cloning’) were nucleofected into 280-NPCs with a piggyBac-transposase-expressing plasmid (see ‘Nucleofection’ for details). Four days later, EGFP⁺ cells expressing the KD constructs were collected by FACS (BD FACSymphony S6). The cells were pelleted and flash-frozen in liquid nitrogen. Protein extraction and auto-western blot were performed by RayBiotech using antibodies against vinculin (R&D systems, MAB6896, clone 728526; diluted 1:25) and WAC (Abcam, AB109486; diluted 1:50).

Gene sets

1. (1)

   SFARI genes: all genes included in the SFARI gene database¹⁴ (3 April 2025 release).

2. (2)

   HC-SFARI genes: SFARI genes with a score of 1.

3. (3)

   Intellectual disability genes: genes from the SysNDD database⁸¹ (accessed 23 August 2025) that were: (i) associated with an intellectual disability term (HP:0001249, HP:0001256, HP:0002342, HP:0010864, HP:0002187); (ii) in the ‘definitive’ category; (iii) not SFARI genes.

Enrichment analyses

Enrichment analyses were performed using one-sided Fisher’s exact tests (fisher.test, base R v.4.3.1), and P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing (p.adjust, method = “BH”, base R v.4.3.1). In bar plots, main bars indicate the OR from a one-sided Fisher’s exact test; error bars indicate the lower bound of the 95% confidence interval on the OR (no upper bound is defined for a one-sided test).

KD of fluorophores

Fluorophore-expressing constructs (‘Molecular cloning’) were nucleofected into 280-iPS cells with a piggyBac-transposase-expressing plasmid (see ‘Nucleofection’ for details). iPS cells were expanded, and cells expressing all three fluorophores were isolated by FACS using a Cytek Aurora CS spectral cell sorter. The sorted iPS cells were further expanded and induced into NPCs (‘Cell culture’). LiOn-based CasRx KD constructs were then nucleofected into these NPCs with a piggyBac-transposase-expressing plasmid. Four days later, cells expressing all four fluorophores (with EGFP indicating expression of the KD constructs) were analysed for fluorescence intensity of the reporters using a Cytek Aurora CS spectral cell sorter. Ordinary two-way ANOVA with Bonferroni’s multiple comparisons test was used to test for significant differences in the geometric mean fluorescence intensity (log₂) across samples and reporters.

KD by shRNA

For each target, pSicoR-based constructs with three different shRNAs were pooled and then nucleofected into 280-NPCs (see ‘Nucleofection’ for details, with modifications as described here). As these constructs are not transposable, the plasmid expressing piggyBac transposase⁶⁹ was omitted, and the amount of pSicoR-based plasmid was increased to a total of 15 μg. The day after nucleofection, EGFP⁺ cells expressing the constructs were isolated by FACS (BD FACSymphony S6) and replated to allow recovery overnight. The following day, RNA was collected using Zymo DNA/RNA Shield. Finally, 3′ RNA-seq was performed by Plasmidsaurus (DGE RNA-Seq service).

DEG analysis of bulk RNA-seq data from shRNA KD NPCs

We performed differential expression analysis using the same analytical pipeline as for our CasRx KD experiments: DESeq2 (v.1.42.0) was used with the model ~targetGene to calculate differential expression relative to NTC shRNAs (designed to target non-expressed transcripts luciferase and lacZ). For transcriptome-wide analyses, P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing for the set of 14,601 genes with baseMean > 10 in NTC cells. DEGs were defined as genes that were differentially expressed with adjusted P < 0.05.

Analysis of gene expression data from ref. ²⁰

Gene expression counts from ref. ²⁰ were downloaded from GEO (GSE107878), and differential expression analysis was performed with the same analytical pipeline as for our CasRx KD experiments: DESeq2 (v.1.42.0) was used with the model expression ~ targeted gene. DEGs were identified relative to control cells for both neurons and iPS cells. Only genes with mean expression greater than 10 in control samples were considered to be potential DEGs, and P values were adjusted accordingly using Benjamini–Hochberg correction. DEGs were defined as genes that were differentially expressed with adjusted P < 0.05.

Cross-enrichment analysis

We quantified the overlap between DEGs identified from: (1) our CasRx KD experiments in NPCs, (2) our shRNA KD experiments in NPCs, (3) knockout iPS cells from ref. ²⁰ and (4) knockout ExNs from ref. ²⁰. For each pairwise comparison of DEGs, enrichment analyses were performed using one-sided Fisher’s exact tests with a universe of genes that were expressed by both sample types. P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing.

Analysis of co-expression with SFARI genes

To determine whether a gene of interest was significantly co-expressed with SFARI genes in our NPC RNA-seq data, we calculated the median Pearson correlation between the expression of the gene of interest and that of each SFARI gene. To assess whether the degree of correlation was significant, we created an empirical null distribution by calculating the median correlation between the gene of interest and 1,000 randomly selected gene sets of the same size as the SFARI gene set, using the universe of NPC-expressed genes (baseMean > 10 in NTCs). A one-sided empirical P value was then calculated as the proportion of null values greater than or equal to the observed median correlation. Regulators with P < 0.05 were considered to be significantly correlated with SFARI genes.

Identification of consistently downregulated DEGs

Consistently downregulated DEGs were defined as genes that exhibited: (1) statistically significant downregulation by at least three of the KDs, and (2) at least 80% concordance in the direction of significant differential expression across all KDs.

Constraint analysis

LOEUF scores (oe_lof_upper) were obtained from gnomad v.2 (ref. ⁸²) for evaluation of genic constraint. For each gene, we calculated a per-gene LOEUF score by calculating the median of the LOEUF values of each gene’s associated transcripts. Although LOEUF is a continuous metric, we used LOEUF < 0.35 as recommended to identify highly constrained genes.

Further support for candidate ASD risk genes

Candidate ASD risk genes (Supplementary Table 8) that were identified in any of the following major sequencing studies are included in Supplementary Table 9: refs. ^{2,52,83,84,85,86}.

Further references provide additional support for some of these candidate ASD risk genes: refs. ^87,88,89,90.

Nucleofection

Nucleofections were performed using a Neon Transfection System (Thermo Fisher Scientific). For NPCs, 2 million cells were nucleofected with 10 μg of transposable plasmid and 2 μg of a plasmid expressing piggyBac transposase⁶⁹ using 100-μl tips. For iPS cells, 1 million cells were nucleofected with 5 μg of transposable plasmid and 1 μg of a plasmid expressing piggyBac transposase⁶⁹ using 100-μl tips. For small-scale nucleofections, 10-μl tips were used, and cell numbers and plasmid amounts were reduced proportionally. Nucleofection of NPCs was performed using the following settings: voltage = 1,600 V, width = 10, pulses = 3. Nucleofection of iPS cells was performed with voltage = 950 V, width = 30, pulses = 2.

NPC proliferation assays

LiOn-based KD constructs (‘Molecular cloning’) were nucleofected into 280-NPCs with a piggyBac-transposase-expressing plasmid (see ‘Nucleofection’ for details). Four days later, the cells were treated for 20 min with dye from a CYTOTRACK Red 628/643 Cell Proliferation Assay Kit (Bio-Rad, 1351205) diluted 1:500 in medium. The cells were then rinsed with phosphate-buffered saline (PBS; Thermo Fisher Scientific, 10010031) and passaged, with one-fifth of each sample plated into a well of each of four separate plates to be used for the subsequent four time points. The remaining cells were analysed by flow cytometry (BD FACSDiva software v.8.0.3) for assessment of the initial CytoTrack intensity on day 0. For the next 4 days, one well of each sample was dissociated each day and analysed by flow cytometry. The CytoTrack intensity values (negative log₂) from the EGFP⁺ cells expressing the KD construct were normalized to those of the EGFP⁻ cells from the same well to control for well-to-well variability in proliferation rates. They were further normalized to their initial CytoTrack intensity from day 0 to control for variability in dye incorporation. The results for each KD are displayed relative to the mean value from the NTC samples, plotted as the mean ± s.e.m. Statistical analysis was performed using ordinary two-way ANOVA with Dunnett’s multiple comparisons test (GraphPad Prism). Full details of sample numbers and statistical results are available in Supplementary Table 10.

Induction of FLICK recombination

To induce expression of Cre and initiate recombination of the FLICK construct, cells or organoids were treated with 2 μg ml⁻¹ doxycycline (Sigma, D9891-1G) for 24 h.

Assessment of FLICK KD efficiency

To determine KD efficiency using FLICK constructs, we analysed three distinct samples per condition. Control 280-iPS cells were nucleofected with constructs encoding the NTC1 control gRNA array or a gRNA array targeting the gene ENSG00000261840, with a plasmid expressing piggyBac transposase⁶⁹ to genomically integrate the FLICK constructs. For induction of KD, the cells were treated with 2 μg ml⁻¹ doxycycline (Sigma, D9891-1G) for 24 h. Cells that had recombined the FLICK construct were then isolated through FACS (BD FACSymphony S6) for EGFP⁺ cells, which were replated in triplicate and grown for a further 2 days. RNA was then isolated using a Direct-zol RNA Microprep Kit (Zymo Research, R2062). To generate cDNA, we reverse transcribed 500 ng of RNA using SuperScript IV VILO Master Mix (Thermo Fisher Scientific, 11756050). We then performed quantitative PCR (qPCR) by combining cDNA at an RNA-equivalent concentration of 1 ng μl⁻¹ (10 ng per 10-μl reaction), primers at a final concentration of 0.5 μM and PowerUp SYBR Green Master Mix for qPCR (Thermo Fisher Scientific, A25742). The qPCR reactions were run in technical triplicate for each sample on a CFX384 Touch Real-Time PCR Detection System (Bio-Rad, 1855484) using the ‘fast cycling mode’ conditions described in the Master Mix protocol. Expression of the KD target gene (ENSG00000261840) was assessed using the ΔΔCt method, with expression first normalized to that of housekeeping gene HPRT1 for each sample and then displayed relative to the average of the NTC control samples, plotted as mean ± s.e.m. Statistical significance was assessed using a two-tailed t-test (GraphPad Prism).

Processing of cerebral organoids for histological analysis

Cerebral organoids were fixed in 4% paraformaldehyde, diluted in PBS from 16% paraformaldehyde (Electron Microscopy Sciences, 15710), at 4 °C overnight. Following fixation, cerebral organoids were rinsed with PBS and then equilibrated in 30% sucrose (w/v in PBS) at least overnight at 4 °C. For large organoids, the 30% sucrose solution was replaced at least once following overnight incubation to ensure complete cryoprotection. For cryosectioning, cerebral organoids were frozen in Tissue-Tek O.C.T. Compound (Sakura, 4583) on dry ice. A Leica CM3050 S cryostat was used to obtain 30-μm sections, which were collected on Fisherbrand Superfrost Plus Microscope Slides (Fisher Scientific, 12-550-15) and stored at −80 °C.

For immunohistochemical analysis, we prepared a blocking buffer by combining 0.3 M glycine (Sigma, G7403-1KG), 1% bovine serum albumin (Sigma, A4503-50G), 0.3% Triton X-100 (Sigma, T8787-50ML), 10% normal goat serum (Fisher Scientific, NC9660079) and PBS. Slides were rinsed in PBS to remove O.C.T. and then incubated in blocking buffer for 1 h at room temperature. Blocking buffer was removed and replaced with fresh blocking buffer containing the primary antibodies of interest, and the slides were incubated in a slide moisture chamber (Avantor, 76278-848) at 4 °C overnight. Slides were then rinsed with PBS three times for 5 mins each time. Secondary antibodies and nuclear stain DAPI (Thermo Fisher Scientific, 62248) were diluted in blocking buffer (1:500 and 1:1,000, respectively) and added to the slides for 30 min at room temperature; the slides were protected from light from this point forward. The slides were then rinsed in PBS three times for 5 min each time. Cover slips (Epredia, 152455) were mounted using Fluoromount-G (Fisher Scientific, OB100-01), which was allowed to dry at least overnight at room temperature. After drying, slides were stored at 4 °C.

Antibodies for immunohistochemistry and immunocytochemistry

The following antibodies were used: PAX6 (BD Pharmingen, 561462, clone O18-1330), diluted 1:300; Nestin (Millipore Sigma, HPA026111-25UL), diluted 1:500; TBR2 (Millipore Sigma, HPA028896-100UL), diluted 1:300; CTIP2 (Abcam, ab18465, clone 25B6), diluted 1:300; EGFP (Aves, GFP-1010), diluted 1:1,000; Goat anti-chicken IgY, Alexa Fluor 488 (Thermo Fisher, A-11039), diluted 1:500; Goat anti-rabbit IgG, Alexa Fluor 488 (Thermo Fisher, A-11008), diluted 1:500; Goat anti-rat IgG, Alexa Fluor 488 (Thermo Fisher, A-11006), diluted 1:500; Goat anti-mouse IgG, Alexa Fluor 555 (Thermo Fisher, A-21424), diluted 1:500; Goat anti-rat IgG, Alexa Fluor 647 (Thermo Fisher, A-21247), diluted 1:500; Goat anti-rabbit IgG, Alexa Fluor 647 (Thermo Fisher, A-21245), diluted 1:500; Goat anti-rabbit IgG, Alexa Fluor 750 (Thermo Fisher, A-21039), diluted 1:500.

Confocal imaging and processing

Cerebral organoid sections were imaged using a Zeiss LSM 980 confocal microscope. Standard image processing was performed using Adobe Photoshop.

Cerebral organoid scRNA-seq

For scRNA-seq, cells were collected from:

• 1 month: 10 organoids, pooled together and collected as a single batch;

• 2 months: 14 organoids collected in 2 batches, with 7 organoids collected together on day 60 and a further 7 organoids collected together on day 61;

• acute KD: 8 organoids, pooled together and collected as a single batch.

For collection of cells, mosaic cerebral organoids were cut into pieces using spring scissors, rinsed twice with PBS to remove debris, and then dissociated using a Neural Tissue Dissociation Kit – Postnatal Neurons (Miltenyi Biotec, 130-094-802) with the enzyme mixes described for brain tissue. Gentle dissociation was performed by trituration with pipette tips of decreasing diameter, starting with wide-orifice tips (Rainin, 30389218) and progressing to standard tips. Cells with recombined FLICK constructs were isolated through FACS (BD FACSymphony S6) for EGFP⁺ cells. Transcriptome-wide scRNA-seq libraries were prepared using Chromium Next GEM Single Cell 3′ Reagent Kits v.3.1 (Dual Index) from 10x Genomics, following user guide CG000316 Rev. A.

Construction of FLICK barcode sequencing libraries

Full details of protocols are available from protocols.io (https://doi.org/10.17504/protocols.io.5qpvoj5pzg4o/v1). Briefly, for quantification of FLICK constructs from the initial pool of iPS cells used to generate mosaic organoids (day 0), RNA was extracted from triplicate samples of 42,000 cells using a Direct-zol RNA Microprep Kit (Zymo Research, R2062). For generation of cDNA, 50–250 ng of RNA was reverse transcribed using a SuperScript IV First-Strand Synthesis System (Invitrogen, 18091050) with an oligo(dT)VN primer containing UMIs. Sequencing libraries were then generated using Q5 Hot Start High-Fidelity 2X Master Mix (NEB, M0494L) and custom primers.

To generate FLICK barcode libraries from scRNA-seq libraries, we used 5 μl of amplified and cleaned cDNA (produced by following the 10x User Guide CG000316 Rev. A to the end of step 2). Barcode-containing cDNA was amplified for sequencing using Q5 Hot Start High-Fidelity 2X Master Mix (NEB, M0494L) and custom primers. Libraries were cleaned using SPRIselect DNA Size Selection Reagent (Beckman Coulter, B23317) and assessed using TapeStation High Sensitivity D1000 reagents (Agilent, 5067-5603 and 5067-5584).

Alignment and quantification of organoid scRNA-seq data

We aligned the scRNA-seq data to the prebuilt Cell Ranger (v.9.0.0) hg38 reference using cellranger-count in feature-barcode mode to perform simultaneous gene expression and guide abundance quantification. Introns were not included in gene expression quantification (–include-introns=false) to avoid quantification of partially degraded transcripts that had been targeted by CRISPR–CasRx. CRISPR guide sequencing library FASTQs for 2-month samples were downsampled to 1.25% of reads using seqtk (v.1.4) to match CRISPR library sequencing depth for 1-month samples and to minimize background for guide assignment using the cellranger-CRISPR guide quantification pipeline. All other parameters were set to default values. Cell Ranger guide assignments were used to filter to only cells with exactly one assigned guide.

Preprocessing and quality control of organoid scRNA-seq data

Long-term and acute KD samples were processed separately but equivalently. CellBender (v.0.3.2) was used for detection and filtering of ambient RNA contamination, whereas scrublet (v.0.2.3) was used for automatic doublet detection and filtering. All other downstream analysis was performed with Seurat (v.5.0.0). High-quality cells were retained according to the following quality control metrics: (1) percent mitochondrial (percent.mt) gene expression < 20% and (2) 500 < number of expressed genes (nFeature_RNA) < 6,000.

Cell type annotation of organoid scRNA-seq data

Cell type annotation was first performed on only NTC cells. NTC cells were extracted to create a control-only Seurat object. This object subsequently underwent normalization (NormalizeData), variable feature identification (FindVariableFeatures), scaling (ScaleData; regressed covariates of percent.mt, percent_ribo, nCount_RNA and nFeature_RNA) and dimensionality reduction (RunPCA and RunUMAP). Integration was performed using harmony (v.1.2.0) across organoid batches, and the integrated object was used for unsupervised clustering (FindNeighbors and FindClusters, resolution values of 0.6 and 1.2). A small number of unsupervised clusters (fewer than 1,000 cells) were filtered owing to low cluster-wide quality control metrics (exceptionally high or low nCount_RNA, high percent.mt). All subsequent two-dimensional visualization of this control-only Seurat object was performed using a UMAP representation of the integrated data. Cell type annotation was performed using canonical brain marker genes. A Seurat object was next constructed for all cells, following the processing pipeline (for instance, normalization, scaling, dimensionality reduction and integration) described above. For cell type annotation, label transfer was performed from the control-only Seurat object (TransferData) to the full Seurat object. Expression patterns of canonical brain marker genes were then used to refine and finalize these cell type annotations.

Analysis of relative abundance per KD in cerebral organoids

The relative abundance of cells for each KD target was determined for each time point (day 0, day 30, day 60 and day 61). For day 0, relative abundance was determined from bulk RNA-seq libraries (‘Construction of FLICK barcode sequencing libraries’) using a custom pipeline (‘Code availability’). For days 30, 60 and 61, relative abundance was determined from scRNA-seq data using the CRISPR guide assignments generated by Cell Ranger (‘Alignment and quantification of organoid scRNA-seq data’). For each sample, results were normalized to the starting proportion at day 0 and displayed relative to NTCs from the same time point. Results were plotted as the mean ± s.e.m. Statistical analysis was performed using repeated measures one-way ANOVA with Dunnett’s multiple comparisons test (GraphPad Prism), considering all long-term samples (days 30, 60 and 61).

Analysis of relative cell type proportions in cerebral organoids

The relative ExN/RG proportion was determined for each KD target at day 30, day 60 and day 61. Results were normalized to the NTC sample from the same time point. Statistical analysis was performed using ordinary one-way ANOVA with Dunnett’s multiple comparisons test (GraphPad Prism), considering all long-term samples (days 30, 60 and 61).

Differential expression analysis of cerebral organoid scRNA-seq data

Differential expression analysis was performed between perturbed cells per guide and cells that received NTCs, stratified by main cell types (radial glia, intermediate progenitor cells, ExNs and inhibitory neurons). Seurat FindMarkers was used to perform differential expression with logfc.threshold = 0, min.pct = 0.10, min.cells.group = 3 and all other parameters set to default values. Enrichment analyses of DEGs were performed using one-sided Fisher’s exact tests, and P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing. The universe was defined as all genes expressed in both comparisons of interest.

Visualization of shared DEGs between organoid cell types

We visualized shared DEGs between organoid cell types (Fig. 3c) using R package igraph (v.2.1.4). The Fruchterman–Reingold layout algorithm (layout_with_fr) was used to determine final node placement.

MAGIC analyses

MAGIC⁹¹ (v.1.1) was used to identify transcription factors and cofactors whose targets were enriched in the DEGs from individual perturbations. The strategy of ‘gene body (promoter to the end of the last exon) plus 1Kb flanking sequence either side of the gene body’ was used. For NPC analyses, the universe was set as all genes with baseMean > 10 in NTC samples. For organoid analyses, the universe was set as all genes expressed in at least 10% of cells of that cell type (Seurat min.pct = 0.10; no log₂ fold-change cutoff), as these were the genes assessed for differential expression.

Re-analysis of REST ChIP–seq data

Data was extracted from Table S2 from ref. ⁹². Mouse genes were mapped to human genes using Orthology.eg.db (v.3.20.0), org.Mm.eg.db (v.3.20.0) and org.Hs.eg.db (v.3.20.0). Enrichment analysis was performed using a one-sided Fisher’s exact test with all BrainSpan genes with normalized expression ≥ 1 as the universe (version ‘RNA-Seq Gencode v10 summarized to genes’; https://brainspan.org; accessed 6 January 2024). P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing. Results are displayed as bar plots, with main bars indicating the OR from a one-sided Fisher’s exact test and error bars indicating the lower bound of the 95% confidence interval on the OR.

GRN construction and benchmarking

A GRN was constructed using bioinformatic tool ARACNe-AP⁹³. The input comprised the RNA-seq data from our 3 NTCs and the successful KD samples with at least 50 DEGs. We used curated lists of transcription factors⁹⁴ and epigenetic regulators⁹⁵ and selected genes expressed in NPCs (baseMean > 10 in NTC samples) as the candidate genetic regulators. To create the GRN, we used default parameters and 100 bootstraps that were then aggregated to form the final network. We then used VIPER⁹⁶ (v.1.36.0) to identify differentially active regulons following standard operating procedures. In brief, we identified regulons with significant differential expression between NTCs and KDs using a t-test, with a null distribution defined by sample permutations.

To benchmark the predictive power of the ARACNe-derived GRNs, we generated further ARACNe networks after individually holding out perturbation samples for 5 target genes that were identified as regulators of at least 250 genes (range: 478–1,219 genes) in their second-degree regulons (that is, their direct targets, which constituted their first-degree regulon, as well as the direct targets of those genes). For each held-out regulator, we assessed the overlap between the predicted positively regulated genes within its second-degree regulon and the set of downregulated DEGs resulting from perturbation of the same regulator, using a one-sided Fisher’s exact test. P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing. P values across the five held-out perturbations were combined using Fisher’s method.

PageRank analysis of GRN

To identify central regulators in our GRN, we converted the GRN to a graph using R package igraph (v.2.1.4) and ran the PageRank algorithm⁹⁷ using igraph method page_rank. We then ordered genes by their PageRank score and ran gene set enrichment analysis⁹⁸ using R package fgsea (v.1.32.4), testing the following ontologies: SFARI genes, HC-SFARI genes, intellectual disability genes (exclusive of SFARI genes) and the 50 Molecular Signature Database Hallmark Pathways⁹⁹ (obtained from https://www.gsea-msigdb.org/gsea/msigdb).

Identification of top candidate regulators for driving convergence in autism

VIPER⁹⁶ (v.1.36.0) was used to identify GRN regulators whose regulons were significantly differentially active across the KDs compared with NTC controls. From these, we prioritized 70 regulators whose regulons were significantly enriched for both SFARI genes (one-sided Fisher’s exact test with Benjamini–Hochberg correction for multiple hypothesis testing) and high PageRank scores (gene set enrichment analysis with Benjamini–Hochberg correction for multiple hypothesis testing).

Further GRN analyses

We calculated the mean normalized degree for the indicated set of regulators and their first-degree regulons (Fig. 4d) using R package igraph (v.2.1.4) by running the function mean(degree(graph, normalized = TRUE)). We then created an empirical null distribution of mean normalized degrees for each indicated set of regulators by creating 1 million subnetworks of the same set size of randomly selected GRN regulators. The relative mean normalized degree was calculated as the observed mean normalized degree for a given set of regulators divided by the mean of the values for the corresponding empirical null distribution. A one-sided empirical P value was then calculated as the proportion of null values greater than or equal to the observed mean normalized degree, with a detection limit of P = 1 × 10⁻⁶.

To evaluate known or predicted protein–protein interactions between the top 70 regulators, we used the STRING database GUI¹⁰⁰. For the universe, we used the full set of genes that were considered as candidate genetic regulators in the GRN: transcription factors⁹⁴ and epigenetic regulators⁹⁵ that were expressed in NPCs (baseMean > 10 in NTC samples). All other parameters were set to default values. For clustering of our 70 regulators into modules, we performed iterative walktrap clustering using R package igraph (v.2.1.4) by running the function cluster_walktrap with default parameters.

Re-analysis of CHD8 ChIP–seq data

For enrichment analysis of CHD8 targets, we used Supplementary Data 1 from ref. ⁴⁵. Enrichment of the indicated gene sets was tested using one-sided Fisher’s exact tests, with the universe set to NPC-expressed genes (baseMean > 10 in NTCs). P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing. Results were displayed using bar plots, with main bars indicating the OR from a one-sided Fisher’s exact test and error bars indicating the lower bound of the 95% confidence interval on the OR.

Identification of female-biased genes

To identify genes that exhibited female-biased expression in the neurotypical brain, we combined DEG call sets from the following work: Table S4 from ref. ¹⁰¹ (identified as significantly female-biased in either the UCLA cohort, BrainVar cohort or through meta-analysis); ref. ¹⁰² (https://storage.cloud.google.com/adult-gtex/bulk-qtl/v8/sb-eqtl/GTEx_Analysis_v8_sbgenes.tar.gz, identified as significantly female-biased in either BRNCTXA or BRNCTXB); Table S4 from ref. ¹⁰³ (identified as significantly female-biased in the brains of all assessed species (human, mice, rats, dogs, and macaques)); and Table S13 from ref. ⁶⁰ (identified as ubiquitously female-biased).

Analysis of snRNA-seq from autistic male individuals

Processed snRNA-seq data from the human frontal cortex was obtained from ref. ⁶⁴. The sex of individual samples was verified on the basis of expression of XIST and Y-linked genes. Only a small number of female samples (one female control sample, one sample from a female participant with idiopathic autism and one sample from a female participant with 15q-duplication-associated autism) were included in this cohort owing to the known male bias of autism. Thus, all downstream analyses focused on the 17 male samples included in this cohort. Differential expression analysis was performed using Seurat function FindMarkers on a cell-type-specific basis between (1) male controls and samples only from male individuals with idiopathic autism; (2) male controls and samples only from male individuals with 15q-duplication-associated autism; and (3) male controls and all samples from male individuals with autism, either idiopathic or 15q-duplication-associated. To control for donor effects, we ran FindMarkers using MAST¹⁰⁴. Significant genes were defined as those with absolute log₂[fold change] > 0.1 and transcriptome-wide P_adj < 0.05.

Plots of sex-differential expression of ZFX

Data were used from refs. ^{60,101,105,106}.

Re-analysis of ZFX targets from refs. ⁵³ and ⁶⁵

ZFX ChIP–seq targets were obtained from Table S3E of ref. ⁶⁵. ZFX targets were also obtained from Table S10 of ref. ⁵³, by San Roman et al., who previously performed integrative analysis of ZFX ChIP–seq data (ENCODE), and ZFX KD RNA-seq data from ref. ⁶⁵. Enrichment was analysed using one-sided Fisher’s exact tests, and P values were adjusted using Benjamini–Hochberg correction for multiple hypothesis testing. For analysis of common ZFX targets from ref. ⁶⁵, the universe was defined as genes expressed in both of the cell lines used for ChIP–seq for which RNA-seq data were also available (C4-2B and MCF-7). For analysis of ZFX targets from ref. ⁵³, the universe was defined as genes expressed in the corresponding cell type from ref. ⁶⁵. Results are displayed using bar plots, with main bars indicating the OR from a one-sided Fisher’s exact test and error bars indicating the lower bound of the 95% confidence interval on the OR.

Software summary

Software used in this work included R (v.4.3.1) and R packages DESeq2 (v.1.42.0), Seurat (v.5.0.0), igraph (v.2.1.4), scrublet (v.0.2.3), harmony (v.1.2.0), Orthology.eg.db (v.3.20.0), org.Mm.eg.db (v.3.20.0), org.Hs.eg.db (v.3.20.0), VIPER (v.1.36.0) and fgsea (v.1.32.4), as well as Cell Ranger (v.9.0.0), CellBender (v.0.3.2), MAGIC (v.1.1), salmon (v.1.9), seqtk (v.1.4) and BD FACSDiva v.8.0.3.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.