Longitudinal Transcriptomic Analysis of Human Cortical Spheroids Identifies Axonal Dysregulation in the Prenatal Brain as a Mediator of Genetic Risk for Schizophrenia

BACKGROUND
Schizophrenia (SCZ) has a known neurodevelopmental etiology, but limited access to human prenatal brain tissue hampers the investigation of basic disease mechanisms in early brain development. Here, we elucidate the molecular mechanisms contributing to SCZ risk in a disease-relevant model of the prenatal human brain.


METHODS
We generated induced pluripotent stem cell-derived cortical spheroids (hCSs) from a large and genetically stratified sample of 14 SCZ patients and 14 age and sex-matched controls (CTRL). The hCSs were differentiated for 150 days, and comprehensive molecular characterization across four time points was carried out.


RESULTS
The transcriptional and cellular architecture of hCSs closely resembled that of 10-24 post-conception week fetal brain, showing strongest spatial overlap with frontal regions of the cerebral cortex. A total of 3,520 genes were differentially modulated between SCZ and CTRL hCSs across organoid maturation, displaying a significant contribution of genetic loading, an over-representation of risk genes for autism spectrum disorder and SCZ, and strongest enrichment for axonal processes in all hCS stages. The two axon guidance genes SEMA7A and SEMA5A, the first a promoter of synaptic functions and the second a repressor, were down and up-regulated in SCZ hCSs, respectively. This expression pattern was confirmed at the protein level and replicated in a large post-mortem sample.


CONCLUSIONS
Applying a disease-relevant model of the developing fetal brain, we identified consistent dysregulation of axonal genes as an early risk factor for SCZ, providing novel insights into the effects of genetic predisposition on the neurodevelopmental origins of the disorder.

organoid maturation (15)(16)(17) or generated brain organoids from healthy control individuals only (18,19).Moreover, betweendonor variability is a major source of gene expression variation in iPSC models, potentially masking case-control differences (20,21), but this critical issue is seldom adequately addressed.Here, we generated iPSC-derived organoids of the cerebral cortex, termed human cortical spheroids (hCSs), from a large, genetically stratified case-control sample (N = 28; control n = 14, SCZ n = 14) and carried out comprehensive transcriptional profiling across 4 time points of hCS maturation.We demonstrate substantial correspondence in patterns of gene expression and cell type composition between the derived hCSs and developing fetal tissue of the cerebral cortex, and reveal a persistent dysregulation of axonal genes in SCZ.

METHODS AND MATERIALS Sample Characteristics
Skin biopsy donors were recruited through the Norwegian TOP (Thematically Organized Psychosis) study.Recruitment procedures, inclusion and exclusion criteria, and clinical assessments for the TOP study as a whole have been described elsewhere (22,23).Fibroblasts were isolated from 14 control donors and 14 donors with SCZ selected based on clinical information and polygenic load.Cases and controls were matched on age and sex, and all except for 3 donors were Europeans.Skin biopsy specimens from all 28 donors were obtained under sterile conditions.Table 1 summarizes the clinical characteristics of the study participants.

Polygenic Risk Score Calculation
Genotyping and polygenic risk score (PRS) calculation were carried out as previously described (24).Briefly, DNA was extracted from blood samples collected at inclusion in the TOP study, and genotyping was performed on Human Omni Express-24 v.1.1 (Illumina).Standard preimputation quality control was performed using PLINK 1.9 (25).PRSice-2 (26) was then used for PRS calculation based on the latest genomewide association study (GWAS) meta-analysis (27).PRS was calculated by aggregating the number of risk alleles carried at each GWAS-defined risk variant weighted by its effect size (28).A GWAS p-value threshold of .5 was used for risk variant selection.

IPSC Derivation and Characterization
Skin fibroblasts from the 28 donors were grown and reprogrammed as previously described (29).Each iPSC line was subjected to rigorous quality control by phenotyping, regular monitoring of morphology, and pluripotency marker expressions at the Norwegian Core Facility for Human Pluripotent Stem Cell Research Centre.KaryoStat GeneChip array (Thermo Fisher Scientific) was used for karyotyping of iPSCs at passages 15 to 20 for digital visualization of chromosome aberrations.

Cortical Spheroids Generation
The iPSCs were differentiated into hCSs following a previously published protocol with proven high reproducibility (30,31

RNA Extraction and Sequencing
Total RNA was extracted from hCSs using the RNeasy Plus Mini Kit (Qiagen).RNA yield was quantified with a NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific), and RNA integrity was assessed with Bioanalyzer 2100 RNA 6000 Nano Kit (Agilent Technologies).Library preparation and paired-end RNA sequencing were carried out at the Norwegian High-Throughput Sequencing Centre (https://www.sequencing.uio.no/).Libraries were prepared with the TruSeq Stranded mRNA kit (Illumina), which involves poly(A) purification to capture coding as well as several noncoding RNAs.The prepared samples were then sequenced in 2 batches on a NovaSeq S4 platform (Illumina) at an average depth of 50 million reads per sample, using a read length of 150 bp and an insert size of 350 bp.

Data Processing
Quality of raw sequencing reads was assessed with FastQC (Babraham Institute).To pass the initial quality control check, the average Phred score of each base position across all reads had to be at least 30.Reads were further processed by cutting individual low-quality bases and removing adapter and other Illumina-specific sequences with Trimmomatic V0.32 (32) using default parameters.HISAT2 (33) was then used to first build a transcriptome index based on Ensembl annotations and next to map the trimmed reads to the human GRCh38 reference genome.To quantify gene expression levels, mapped reads were summarized at the gene level using featureCounts (34) guided by Ensembl annotations.

Cell Type Deconvolution for In Vivo Spatiotemporal Brain Expression Comparison
Computational estimation of cell type abundances (deconvolution) for spatiotemporal comparisons was performed with CIBERSORTx (35)

Weighted Gene Coexpression Network Analysis
Gene coexpression networks were constructed using the WGCNA package (39).Expression data were transformed using the variance stabilizing transformation method in the DESeq2 package (40).Signed coexpression networks were constructed by first calculating the pairwise correlations between gene expression profiles using the Pearson method.Gene modules were identified by hierarchical clustering, considering only modules with at least 25 genes.To relate the identified modules to external variables (stage of hCS differentiation, estimated cell type proportions, and diagnostic status), Pearson's correlation and the corPvalueStudent function were used to test for significant associations between these variables and each module's module eigengene.To identify the major module-specific driver genes, we made use of the concepts of gene significance and module membership, denoting the extent to which a gene's expression profile correlates with an external trait of interest and the module eigengene of a given module, respectively.Module-specific driver genes were defined as the genes with the largest gene significance plus module membership values.To functionally annotate the identified gene expression modules, Gene Ontology (GO) overrepresentation tests were carried as described next, using only module genes with module membership .0.7 as input.

Differential Expression Analysis and GO Enrichment Tests
Before conducting differential expression (DE) analyses, lowly expressed and nonexpressed genes were filtered out by requiring more than 1 count per million in at least 14 samples (the smallest group sample size).As KaryoStat analysis revealed chromosomal aberrations in some of the iPSC lines, all genes located within affected regions were excluded before analysis (580 genes in total) (Table S1).After applying these filtering steps, a total of 16,225 genes were retained.Time-course DE analyses were conducted using the limma-voom framework (41) adjusting for baseline differences between cases and controls at the iPSC stage, such that only genes displaying variable Prenatal Axonal Dysregulation in Schizophrenia Biological Psychiatry April 1, 2024; 95:687-698 www.sobp.org/journalBiological Psychiatry differences between control and SCZ samples across time were considered.As principal component analysis and hierarchical clustering revealed 3 outlier samples at the iPSC stage (Figure S1), the voomWithQualityWeights function was applied to down-weight the contributions from these samples, thereby reducing noise while preserving statistical power (42).Demographic characteristics such as age, sex, and ethnicity were not significantly different between cases and controls (Table 1); therefore, only sequencing batch was included as a covariate in the linear models.A differentially expressed gene (DEG) was considered significant if the false discovery rate (FDR) was , .1.GO enrichment tests of significant gene sets were conducted with clusterProfiler (43)."Biological processes" was selected as the ontology of interest, and a GO term was considered significantly enriched if the FDR was , .1.

Concordance Analysis of Axonal Genes
The CommonMind Consortium (CMC) dataset ( 44) was used to test whether the identified axonal genes with consistent expression trajectory showed concordant direction of effect in postmortem dorsolateral prefrontal cortex samples from cases and controls.Only samples of European origin provided by the CMC were included (n = 429), comprising 215 control and 214 SCZ donors.Prefiltering and DE analysis were carried out as described earlier, adjusting for age, reported gender, brain bank, and postmortem interval.

Derivation and Molecular Characterization of hCSs
We derived hCSs from well-characterized SCZ and control donors using an established protocol designed to generate reproducible three-dimensional neural structures patterned after the cerebral cortex (30,31).The hCSs were cultured for 150 days, and 4 stages across differentiation (iPSC, D30, D90, and D150) representing key neurodevelopmental transitions were selected for transcriptional profiling (112 samples in total) (Figure 1A, B).The donors were matched on age and sex, and the majority (25 of 28, 89%) were Europeans (Table 1).Importantly, SCZ donors were selected based on PRS for SCZ, and mean PRS was significantly higher in the SCZ group (mean [SD] = 0.33 [0.5] vs. 20.83[1.52]; p = .026)(Figure 1C), which exhibited less variability than the control group, indicating a homogeneous and disease-relevant genetic background of the SCZ donors.The hCS size increased steadily and reached a mean volume of 7.85 mm 3 (2.7 mm in diameter) at D150 (Figure 1D).Expression trajectories of broad markers for cortical maturation followed expected patterns (45), with pluripotency markers peaking first, followed by neuronal markers and then astrocyte markers, which peaked at D150, reflecting the transition from neurogenesis to gliogenesis (Figure 1E).Four markers (SOX2, SOX1, PAX6, and POU3F2) displayed DE between cases and controls (Figure 1E).

Correspondence Between hCSs and Human Primary Brain Tissue
To determine whether the selected stages of hCS differentiation exhibited distinct transcriptional profiles, we performed principal component analysis.The first two principal components (explaining 64% of the variance) showed clear separation between stages, indicating that hCS differentiation is the major contributor to gene expression variability (Figure 2A).The less pronounced separation seen between D90 and D150 suggests that hCSs already reach transcriptional maturity at D90. Stratifying samples by diagnostic status revealed a progressive separation between control and SCZ samples (Figure S1A).Disease-specific differences thus appear to manifest more strongly as hCSs mature.The principal component analysis also revealed 3 outlier samples within the iPSC group (Figure 2A), and this was confirmed by hierarchical clustering (Figure S1B, C).These samples were either removed or otherwise accounted for in subsequent analyses (see Methods and Materials).
To assess the extent to which the derived organoids transcriptionally match primary human brain tissue, spatiotemporal comparisons were conducted using BrainSpan, an anatomically comprehensive atlas of human brain development (37,38).Consistent with previous reports (18,19,46), hCS profiles overlapped with fetal brain tissue at 10-to 24 postconception weeks, beginning at in vitro stage D30 and becoming more prominent by D150 (Figure 2B).The iPSC stage did not map to any temporal period in vivo, likely reflecting the fact that the BrainSpan reference does not include the embryonic phase of prenatal development (before postconception week 8) when brain cells may still retain pluripotency (11) (Figure 2B).Spatial comparisons showed that hCSs gradually assumed region-specific cortical identity, displaying strongest overlap with frontal regions of the cerebral cortex (medial frontal cortex, dorsolateral prefrontal cortex, and orbitofrontal cortex) and weakest overlap with the cerebellar cortex at D150 (Figure 2C).Assessing differences between control and SCZ hCSs, we observed an unexpected temporal pattern in which SCZ hCSs initially matured faster than control organoids, but this trend was reversed by D150 (Figure 2D).Interestingly, this phenomenon of accelerated differentiation has also been observed in NPCs derived from patients with bipolar disorder transitioning toward neuronal identity (47).
To characterize the cellular composition of hCSs compared to primary tissue, computational estimation of cell type proportions (deconvolution) was conducted using human reference signatures for 6 broad and hCS-relevant cell populations (see Methods and Materials).At D30, most iPSCs had differentiated into NPCs and mature neurons, marking an already advanced neurogenesis phase in vitro (Figure 2E, left).The neuronal fraction reached a plateau at D90, accompanied by a corresponding decrease in NPCs and a steady increase in astrocyte abundance, which peaked at D150 (Figure 2E, left).We observed a remarkable degree of correspondence between in vitro and in vivo cell type trajectories within the relevant time periods (Figure 2E, shaded areas), with neurons, astrocytes, and NPCs constituting the main cell types at hCS stage D150 as well as postconception weeks 19 to 24 in vivo, with oligodendrocytes appearing later (Figure 2E).Overall, these findings demonstrate that our organoids preserve key cell type

Gene Coexpression Networks in Maturing hCSs
To further elucidate the transcriptional landscape of maturing hCSs, we performed weighted gene coexpression network analysis (39).We identified 11 gene expression modules, 5 of which correlated significantly with hCS differentiation (p , 2 3 10 29 ) (Figure 3A; Figure S2; Table S2).Modules 2 (93 genes) and 3 (3072 genes) followed a downward trajectory and were primarily enriched for biological processes (GO terms) related to protein translation and cell proliferation, respectively (Figure 3A, B; Table S3).On the other hand, modules 4 (675 genes) and 10 (3009 genes) exhibited an upward trajectory comprising genes mostly related to macroautophagy and synaptic processes, respectively (Figure 3A, B; Table S3).This gradual loss of cell proliferation capacity combined with the reciprocal increase in synaptic signaling ability reflects the progressive differentiation of iPSCs into neurons (Figures 2E  and 3A).Module 8 (315 genes) also displayed an upward trajectory similar to modules 4 and 10, but it was not significantly enriched for any biological process.Module 4 was significantly downregulated in SCZ hCSs at D150, with the antiapoptotic gene BCL2L1 (BCLX) showing the strongest association (Figure 3C; Figure S3; Table S4).This gene is known to promote neuronal cell survival ( 48) and prevent axonal degradation (49).A separate gene coexpression analysis was carried out on the data at D150 alone to investigate the association between gene expression modules and genetic loading.We found that the PRS was associated with a large module enriched for genes involved in neuronal processes and Wnt signaling, suggesting that predisposing genetic factors play an important role in disease-related cellular processes (Figure S4).

Persistent Dysregulation of Axonal Genes in SCZ hCSs
As SCZ might also be strongly influenced by individual genes not fitting neatly into a specific coexpression module, we performed transcriptome-wide time-course DE analysis adjusting for baseline differences at the iPSC stage.We identified 654, 1980, and 2142 DEGs in SCZ hCSs (FDR , .1) at D30, D90, and D150, respectively, with partial overlap between the gene sets (Figure 4A, B; Table S5) and a considerable contribution of PRS loading, particularly during early stages of differentiation (Figure S5).Mapping these DEG sets onto risk genes for 5 major neuropsychiatric, neurodevelopmental, and neurodegenerative disorders, we found the strongest associations for autism spectrum disorder (ASD) (FDR = 3.7 3 10 29 ) and SCZ (FDR = 1.8 3 10 24 ) at D150 (Figure 4C; Table S6; Supplemental Methods).The DEG sets were further annotated based on enrichment for biological processes.The GO terms axonogenesis and/or axon development were the processes with strongest enrichment in all DEG sets (Figure 4D; Table S7), suggesting a consistent axonal

Prenatal Axonal Dysregulation in Schizophrenia
Biological Psychiatry deficit in SCZ throughout hCS maturation.Of the 178 axonal DEGs underlying this enrichment, 21 were consistently upregulated or downregulated in SCZ hCSs from D30 to D150 (Figure S6).To investigate whether the dysregulation of these axonal genes is also present after disease onset, DE analysis was conducted using human primary brain tissue (dorsolateral prefrontal cortex) from 214 control and 215 SCZ donors of European ancestry from the CMC dataset (44).Of the 21 axonal genes with consistent case-control differences throughout hCS maturation, we found that 11 were also significantly dysregulated in CMC (FDR , .1), and 10 of these showed concordance with respect to direction of effect (Figure 4E).

Aberrant Semaphorin Modulation in SCZ hCSs
The genes showing concordant effects in hCSs and postmortem brain tissue included SEMA7A and SEMA5A, both of which encode members of the semaphorin family of proteins known for their pivotal role in axon guidance (50).Importantly, SEMA7A was downregulated in SCZ, while SEMA5A was upregulated (Figure 5A), and this pattern was confirmed in hCSs both by quantitative polymerase chain reaction and at the protein level by immunostaining (Figure 5B-D; Supplemental Methods).As semophorin genes show preferential expression in some cell types (https://www.proteinatlas.org),we investigated whether or not the SCZrelated changes in SEMA7A and SEMA5A expression were driven by differences in cell type composition.We first generated a high-confidence gene expression signature by performing single-cell RNA sequencing on 3 of the control hCS lines differentiated for 180 days (Supplemental Methods).Data integration and annotation resulted in the identification of 11 cell type clusters displaying two primary trajectories, one for neuronal development and one for glial development, both originating from NPCs (Figure 5E).Using this single-cell RNA sequencing signature as prior information, we employed a Bayesian method (51) to reconstruct cell type composition in hCSs and found significant differences between SCZ and control hCSs in the estimated fractions of mature neurons, oligodendrocyte precursor cells and radial glial cells at D150 (Figure S7).Importantly, a strong positive correlation was found between SEMA5A and oligodendrocyte precursor cell fraction, suggesting that the increased expression of SEMA5A in SCZ hCSs is partly due to an increased oligodendrocyte precursor cell proportion (Figure 5F; Figure S7).In addition to their axonal functions, semaphorins have been shown to be involved in immunological processes (52,53).Using enzyme-linked immunosorbent assay, we quantified the extracellular levels of 5 cytokines (interleukin [IL]-1b, IL-10, IL-8, IL-18, and tumor necrosis factor a) known to be modulated by SEMA7A and SEMA5A (54-56), and found no significant differences between SCZ and control hCSs at any stage of maturation (Figure S8), suggesting that the association of these two genes with SCZ is primarily related to their axonal functions.

DISCUSSION
We comprehensively assessed maturing hCSs derived from a large and genetically stratified sample of control and SCZ donors and found that the hCSs recapitulated the transcriptional architecture of human fetal tissue of the frontal cortex and preserved key cellular transitions.We further identified transcriptional deficits in axonal genes as an early neuropathological mechanism in SCZ.
According to the neurodevelopmental paradigm of SCZ, disturbances in early brain development may later contribute to clinical manifestations, when maturing neural systems and circumstantial pressures pose increasing demands on brain functioning that expose the preexisting abnormalities (57,58).The fetal phase is particularly vulnerable.During this period, the repertoire of neurons found in the adult brain is formed, and most anatomical structures are established (11).It is also the molecularly most dynamic period of human brain development (59), making it sensitive to pathological genetic influences, as underscored by the observation that multiple genes within SCZ-associated risk loci are preferentially expressed in the prenatal brain (60)(61)(62)(63)(64).The importance of the prenatal period is further attested to by the high degree of genetic overlap between SCZ and ASD, especially for rare coding (65,66) and copy number variants (67,68), supporting the idea that SCZ lies on a neurodevelopmental continuum with ASD and other disorders of the developing brain (2).Consistent with this, we found that the DEGs associated with SCZ across hCS maturation were enriched for ASD and SCZ risk genes.The stronger association with ASD likely reflects the fact that the gene reference used for ASD included all types of genetic variation, while the SCZ gene reference included genetic evidence only from common variants.
Our results point to several axonal genes as potential early contributors to SCZ risk, including SEMA7A and SEMA5A, which encode two members of the semaphorin family of proteins originally discovered as guidance cues for developing axons (69,70).SEMA5A is a known ASD susceptibility gene (71,72), and SEMA3A, another member of the family, has been implicated in SCZ (73,74).Importantly, SEMA7A has been found to promote axon growth and neurogenesis (75,76), while SEMA5A has been identified as a negative regulator of synapse formation (77).In both our hCSs and the postmortem CMC sample, SEMA7A was downregulated in SCZ, while SEMA5A was upregulated, suggesting that dysregulation of these two genes may suppress axonal processes in SCZ through distinct but mutually reinforcing mechanisms.
Burgeoning evidence has demonstrated that semaphorins also exert immunomodulatory effects in the nervous system (54), but we found no differences between control and SCZ hCSs in secreted levels of key cytokines.While this suggests that our SEMA findings are related to their axonal roles, an immunological contribution cannot be ruled out due to the absence of peripheral immune cells in our hCS cultures.Overall, our results are in line with a similar study showing downregulation of pathways related to neurodevelopment and synaptic biology in SCZ (15) and with the identification of genes implicated in synaptic biology in the most recent SCZ GWAS (10).
A well-known challenge with iPSC models is that a high degree of gene expression variability between individuals may obscure disease-specific variation and confound reproducible disease modeling (21,78).A strength of our study was the use of a large number of donors (N = 28) with genetically homogeneous SCZ cases as measured by a high group-level PRS, which captures an individual's overall genetic loading for SCZ (79).This strategy is an effective way to reduce the negative influence of donor heterogeneity and increase statistical power (20,80), suggesting that our findings reveal genuine diseaserelevant mechanisms and not spurious genetic effects.However, it should be noted that although the derived hCSs mimic the tissue architecture and transcriptional trajectories of developing fetal cortex, they lack important features that are expected to affect disease progression, including vascularization and physiological and environmental inputs (6,12).Thus, our hCS model primarily captures abnormal disease processes stemming from predisposing genetic factors.

Conclusions
Using reproducible organoids of the frontal cortex, we modeled the neurodevelopmental mechanisms of SCZ in the critical period of prenatal brain development and identified persistent transcriptional aberrations in axonal genes.Our results point to novel disease mechanisms and further advance understanding of the neurodevelopmental origins of SCZ.

PrenatalFigure 1 .
Figure 1.Derivation and characterization of hCSs.(A) Experimental design of the study.(B) Main growth factors and media used for neural induction and differentiation.DM and SB are SMAD inhibitors; XAV is a WNT inhibitor.(C) PRSs in CTRL (n = 12) and SCZ (n = 10) donors.The 3 non-European donors were excluded as PRS estimates from current genome-wide association studies are more accurate in European populations.PRS data were not available for 2 CTRL donors and 1 SCZ donor.*p , .05,two-sided Welch test.(D) Mean volume of all hCSs (n = 28) at each stage of maturation.Error bars indicate SD. (E) Transcriptional trajectories of cell type-specific marker genes showing expected gene expression patterns.*p , .05,**p , .01,two-sided Mann-Whitney test.BDNF, brain-derived neurotrophic factor; CTRL, control; EGF, epidermal growth factor; FGF-2, fibroblast growth factor 2; hCS, human cortical spheroid; iPSC, induced pluripotent stem cell; NT3, neurotrophin-3; PRS, polygenic risk score; SCZ, schizophrenia.

Figure 2 .
Figure 2. Transcriptional overlap between hCSs and human in vivo brain tissue from BrainSpan.(A) PC analysis plot showing good sample clustering by hCS differentiation stage.Note the 3 outlier samples within the iPSC group.(B) Gene expression overlap between hCS profiles and in vivo temporal profiles showing strongest overlap with in vivo periods PCW 10-24 (shaded area).(C) Gene expression overlap between hCS profiles and in vivo spatial profiles.Strongest overlap was seen with the in vivo cortical regions MFC, DFC, and OFC.Weakest overlap was seen with CBC.(D) Gene expression overlap between hCS profiles and in vivo temporal profiles stratified by diagnostic status.Significant differences, and in the opposite direction, seen for hCS stages D30 and D150.*p , .05,two-sided Welch test.(E) Computational deconvolution showing strong concordance in estimated cell type fractions across in vitro and in vivo stages with strongest overlap (shaded area).Note the decreasing pattern of the neuronal fraction in vivo after PCW 16-18, which is not due to a decrease in absolute cell numbers, but rather a relative decrease in the neuronal proportion as the astrocyte and oligodendrocyte fractions increase.A1C, primary auditory cortex; AMY, amygdala; CBC, cerebellar cortex; CTRL, control; DFC, dorsolateral prefrontal cortex; hCS, human cortical spheroid; HIP, hippocampus; IPC, inferior parietal cortex; iPSC, induced pluripotent stem cell; ITC, inferolateral temporal cortex; M1C, primary motor cortex; MD, mediodorsal nucleus of thalamus; MFC, medial frontal cortex; NPC, neural progenitor cell; OFC, orbitofrontal cortex; PC, principal component; PCW, postconception week; S1C, primary somatosensory cortex; SCZ, schizophrenia; STC, superior temporal cortex; STR, striatum; V1C, primary visual cortex; VFC, ventrolateral prefrontal cortex.

Figure 3 .
Figure 3. Coexpression networks in differentiating human cortical spheroids.(A) Expression of top 5 significant The module eigengene is the representative of the gene expression profiles in the module.(B) Functional annotation of top 5 modules by Gene Ontology term enrichment (biological processes).The size of the bubbles corresponds to the 2log10 (p value) of each Gene Ontology term.M8 was not significantly enriched for any biological process.(C) Module trajectories stratified by diagnostic status showing a significant difference between SCZ and CTRL for M4 at human cortical spheroid stage D150.*p , .1, two-sided Mann-Whitney test.CTRL, control; ER, endoplasmic reticulum; iPSC, induced pluripotent stem cell; M, module; SCZ, schizophrenia.

Figure 4 .
Figure 4. Consistent dysregulation of axonal genes in SCZ hCSs and postmortem brain tissue.(A) Bar plot showing upregulated and downregulated DEGs in SCZ across hCS maturation.(B) Overlap between DEG sets for each hCS stage.(C) Enrichment of identified SCZ DEGs in risk genes for 5 neuropsychiatric, neurodevelopmental, and neurodegenerative disorders, showing strongest correlation for ASD and SCZ.Numbers in each box indicate strength of enrichment.The disease gene sets were obtained from either the latest and largest GWAS for each disorder or from comprehensive databases as indicated in parentheses.(D) Gene Ontology enrichment analysis (biological processes) of identified DEGs for each hCS stage.Bubble size corresponds to 2log10 (p value) of enriched Gene Ontology terms.(E) Concordance of consistent axonal genes in hCSs and postmortem brain tissue (dorsolateral prefrontal cortex) using the CMC dataset (SCZ n = 214, CTRL n = 215).Purple bars indicate genes with concordant direction of effect.AD, Alzheimer's disease; ASD, autism spectrum disorder; BD, bipolar disorder; CMC, CommonMind Consortium; CTRL, control; DEG, differentially expressed gene; FDR, false discovery rate; GWAS, genome-wide association study; hCS, human cortical spheroid; iPSC, induced pluripotent stem cell; PD, Parkinson's disease; SCZ, schizophrenia.

Figure 5 .
Figure 5. Aberrant semaphorin modulation in SCZ hCSs.(A) RNA expression of SEMA7A and SEMA5A showing concordant dysregulation in hCSs and postmortem brain tissue measured by RNA sequencing.(B) Confirmation by quantitative polymerase chain reaction of dysregulated SEMA7A and SEMA5A expression in SCZ hCSs.(C, D) Immunohistochemistry of SEMA7A and SEMA5A in representative D150 hCSs (C) and quantification of protein expression levels (D).(E) UMAP plot showing the 11 cell type clusters identified by single-cell RNA sequencing of 3 CTRL hCSs.Numbers correspond to the numbered labels in panel (F).The left arm represents the neuronal trajectory, while the right arm represents the glial trajectory.(F) Single cell-based estimation of cell type fractions and correlation with SEMA7A and SEMA5A expression at D150. *p , .05,**p , .01,two-sided Mann-Whitney test.a.u., arbitrary unit; CMC, CommonMind Consortium; CTRL, control; hCS, human cortical spheroid; iPSC, induced pluripotent stem cell; mRNA, messenger RNA; NPC, neural progenitor cell; OPC, oligodendrocyte precursor cell; RG, radial glial (cell); SCZ, schizophrenia.