Human forebrain organoids-based multi-omics analyses reveal PCCB’s regulation on GABAergic system contributing to schizophrenia

Identifying genes whose expression is associated with schizophrenia (SCZ) risk by transcriptome-wide association studies (TWAS) facilitates downstream experimental studies. Here, we integrated multiple published datasets of TWAS (including FUSION, PrediXcan, summary-data-based Mendelian randomization (SMR), joint-tissue imputation approach with Mendelian randomization (MR-JTI)), gene coexpression, and differential gene expression analysis to prioritize SCZ candidate genes for functional study. Convergent evidence prioritized Propionyl-CoA Carboxylase Subunit Beta (PCCB), a nuclear-encoded mitochondrial gene, as an SCZ risk gene. However, the PCCB’s contribution to SCZ risk has not been investigated before. Using dual luciferase reporter assay, we identified that SCZ-associated SNP rs35874192, an eQTL SNP for PCCB, showed differential allelic effects on transcriptional activities. PCCB knockdown in human forebrain organoids (hFOs) followed by RNA-seq revealed dysregulation of genes enriched with multiple neuronal functions including gamma-aminobutyric acid (GABA)-ergic synapse, as well as genes dysregulated in postmortem brains of SCZ patients or in cerebral organoids derived from SCZ patients. The metabolomic and mitochondrial function analyses confirmed the deceased GABA levels resulted from reduced tricarboxylic acid cycle in PCCB knockdown hFOs. Multielectrode array recording analysis showed that PCCB knockdown in hFOs resulted into SCZ-related phenotypes including hyper-neuroactivities and decreased synchronization of neural network. In summary, this study utilized hFOs-based multi-omics data and revealed that PCCB downregulation may contribute to SCZ risk through regulating GABAergic system, highlighting the mitochondrial function in SCZ.


Introduction
Schizophrenia (SCZ) is a complex polygenic psychiatric disorder with risk contributed by environmental and genetic factors 1 . Genetic studies such as genome-wide association studies (GWAS) have identi ed hundreds of common single nucleotide polymorphisms (SNPs) associated with SCZ 2,3 . Most of the SCZ-associated SNPs are non-coding variants located in regulatory DNA elements [4][5][6] , suggesting that gene expression mediates the connection between genetic variants and SCZ phenotypes 7 . Identifying genes whose expression is associated with SCZ phenotypes facilitates discovering SCZ risk genes for downstream functional studies.
By integrating SCZ GWAS and brain expression quantitative trait loci (eQTL) data, several approaches which are collectively described as transcriptome-wide association studies (TWAS) have been used to identify SCZ risk genes. These TWAS approaches, including FUSION 8-10 , PrediXcan 11 , summary-data-based Mendelian randomization (SMR) 2,10,12,13 , and joint-tissue imputation approach with Mendelian randomization (MR-JTI) 14 , aimed to identify the association between predicted gene expression and SCZ risk. Though MR-JTI could improve gene expression prediction performance in TWAS and provide a causal inference framework 15 , experimental validation is still needed.
Here we integrated results from MR-JTI 14 and other published SCZ TWAS datasets 2,8−13 to prioritize SCZ risk genes for functional study.
Through these procedures, we identi ed Propionyl-CoA Carboxylase Subunit Beta (PCCB), a protein-coding gene that plays important roles in mitochondrial metabolism 16,17 , as an SCZ risk gene with the most supporting evidence in our analysis. However, the PCCB's contribution to SCZ risk has not been investigated before. Using human forebrain organoids (hFOs), three-dimensional cell cultures that capture key aspects of human brain 18 , we found that PCCB knockdown in hFOs resulted into SCZ pathology-related cellular phenotypes. We also identi ed that SCZ-associated SNP rs35874192 may regulate PCCB expression, supporting that SCZ-associated common genetic variants may regulate PCCB expression which mediates the genetic effects on SCZ risk.

Results
PCCB is prioritized as a promising SCZ risk gene To obtain reliable SCZ risk genes for downstream functional study, we integrated multiple published TWAS datasets (Table S1) to prioritize genes with su cient supporting evidence. We also checked whether the prioritized genes are located in SCZ risk-associated gene coexpression module or dysregulated in postmortem SCZ brains. These analyses prioritized PCCB, GATAD2A, and GNL3 as the top three SCZ risk genes (Table 1). Notably, PCCB was also identi ed as an SCZ risk gene in the gene-based MAGMA analysis 19 . Moreover, PCCB is located in the gene coexpression module (M2) that is downregulated in SCZ based on the PsychENCODE data 10 . PCCB was also found to be nominally downregulated in postmortem SCZ brains (P = 0.01, FDR = 0.14) in the CommonMind data 20 . These lines of evidence suggested that PCCB expression mediated the genetic effects on SCZ risk. Therefore, we focused on studying how PCCB contributes to SCZ risk in this study. Since PCCB expression is genetically associated with SCZ, we investigated the functional impacts of SCZ-associated SNPs on PCCB expression. Based on the TWAS results used in this study, we retrieved the top SNPs (rs7432375, rs7427564, rs527888, rs66691851) and their linkage disequilibrium (LD) SNPs that were associated with PCCB expression. To narrow down to the putatively causal variants, we focused on those eQTL SNPs (eSNPs) that are likely to affect PCCB expression in the brain. Since opening chromatin facilitates gene expression activation, we used brain ATAC-seq data from the PsychENCODE consortium 21 to identify eSNPs located in active transcription regions. By integrating PsychENCODE ATAC-seq data and SNP annotation information from the Roadmap Epigenetics Consortium 22 , we prioritized three eSNPs (rs35874192, rs900818, rs7349597) ( Table 2) that are located in genomic regions strongly suggested as enhancers or promoters in the human brain tissues or neural cell cultures. We then performed dual luciferase reporter assay (DLRA) in both hNPCs and SH-SY5Y cell lines to validate the regulatory effects of the three eSNP-containing DNA elements. For each eSNP, 50 base pairs (bp) eSNP-containing DNA fragment was synthesized and cloned into the upstream of PCCB promoter in the PGL3-basic luciferase reporter vector (Fig. S1). In DLRA, the eSNP rs35874192 (G/C) showed allelic effects on transcriptional activities in both hNPCs and SH-SY5Y cell lines, with SCZ-associated allele C corresponding to a lower gene expression (Fig. 1A). Notably, the directions of allelic effects of rs35874192 were consistent with eQTL patterns detected in brain tissues from the GTEx 23 and BrainSeq consortium 24 (Fig. 1A). While the eSNPs rs900818 and rs7349597 had no differential allelic effects on transcriptional activities (Fig. 1B, C).

PCCB knockdown in hFOs affects expression of genes enriched in GABAergic synapse
According to the TWAS and DLRA results, lower PCCB expression is associated with increased SCZ risk. We established PCCB knockdown and control human induced pluripotent stem cells (hiPSCs, U2F) using CRISPR interference (CRISPRi). In CRISPRi, one guide RNA (gRNA) sequence targeting PCCB (PCCB-G1) and one non-targeting control gRNA were designed. The established PCCB knockdown and control hiPSCs were then used to generate hFOs ( Fig. 2A, B) to investigate the functional impacts of PCCB knockdown. On day 60 of organoid culture ( Fig. 2A, C, D), PCCB knockdown and control hFOs were used for RNA-sequencing (RNA-seq). Differential gene expression analysis identi ed 2326 differentially expressed genes (DEGs) [false discovery rate (FDR) < 0.05] between the PCCB knockdown and control hFOs (Table S2). Among the 2326 DEGs, 1099 genes were upregulated and 1227 genes were downregulated in PCCB knockdown hFOs (Fig. 2E, F).

PCCB -induced DEGs in hFOs are enriched with SCZ-related genes
To explore PCCB's connection to SCZ, we evaluated the enrichment of 1079 PCCB-induced DEGs in hFOs with SCZ-related gene sets. The rst SCZ gene set was 4096 differentially expressed protein-coding genes between postmortem brains of 559 SCZ patients and 936 controls from the PsychENCODE consortium 10 . The second SCZ gene set was 2809 DEGs between cerebral organoids (6 months) derived from eight SCZ patients and eight controls from the Kathuria et al. study 26 . We found that the 1079 PCCB-induced DEGs were signi cantly overlapped with genes dysregulated in PsychENCODE SCZ brains (Overlapped genes = 282, P = 5.84E-04) and SCZ patient-derived cerebral organoids (Overlapped genes = 255, P = 1.92E-14) (Fig. S3A). We also found that the 1079 PCCB-induced DEGs were signi cantly (P adjust =7.09E-3) overlapped with genes reported in SCZ GWAS from the FUMA analysis 27 (Fig. S3B). These results suggested that PCCB may contribute to SCZ risk by affecting expression of genes related to SCZ (Table S2).
PCCB knockdown in hFOs decreases GABA level by reducing tricarboxylic acid cycle and leads to mitochondrial dysfunction PCCB encodes the β subunit of the propionyl-CoA carboxylase, a mitochondrial enzyme involved in the catabolism of propionyl-CoA 28 . PCCB mutation has been reported to impair mitochondrial energy metabolism by disrupting the tricarboxylic acid (TCA) cycle 16 . We would expect mitochondrial dysfunction caused by PCCB knockdown. Indeed, we found several mitochondrial genes that function in cellular oxidative phosphorylation, including MT-ND2, MT-ND5, and MT-CYB, were downregulated in both PCCB-G1 and PCCB-G2 hFOs (Table S2), which was further validated by the RT-qPCR analysis (Fig. 4A). Adenosine triphosphate (ATP) and reactive oxygen species (ROS) detection assays showed that PCCB knockdown reduced ATP generation and increased ROS levels in hFOs (Fig. 4B, C), indicating the mitochondrial dysfunction caused by PCCB knockdown.
Since GABA metabolism involves a route from α-ketoglutarate (α-KG) generated by the TCA cycle to succinate via glutamate, GABA, and succinic semialdehyde 29,30 , we examined whether PCCB knockdown decreased GABA level by inhibiting the TCA cycle. We performed enzyme linked immunosorbent assay (ELISA) and con rmed that PCCB knockdown decreased succinyl-CoA (SCOA) and α-KG (Fig. 4D, E), two key metabolites that connect the GABA shunt and TCA cycle 29,30 . Considering that α-KG is an upstream metabolite that could be converted to GABA, we added α-KG (10 µg/ml) into culture media of hFOs and found a restored GABA level (Fig. 4F) in PCCB knockdown hFOs. These results indicated that PCCB knockdown decreased GABA level by reducing TCA cycle (Fig. 4G).
PCCB knockdown in hFOs leads to abnormal electrophysiological activities Since GABA, the major inhibitory neurotransmitter in the brain 31 , was decreased in PCCB knockdown hFOs, we used multielectrode array (MEA) recording assay to test whether PCCB knockdown affected neuroactivities. The PCCB knockdown and control hFOs at day 160 were seeded on Matrigel-coated 24-well MEA plate (Fig. 5A). After 7 days of culture, electroactivities of hFOs were recorded (Fig. 5B, C, D). We found that PCCB knockdown in hFOs led to increased number of spikes (Fig. 5E) and mean neuron ring rate (Fig. 5F), suggesting a hyper neuroactivity after PCCB knockdown in hFOs. However, PCCB knockdown decreased the synchronization of the neural network (Fig. 5G).
Since hyper neuroactivity and decreased synchronization of neural network in SCZ brains have been reported by the electroencephalography and magnetoencephalography 32 , these results supported that PCCB knockdown led to abnormal electrophysiological activities that link to SCZ phenotypes.

Discussion
SCZ is a polygenic psychiatric disorder with risk contributed by multiple genes. Identifying genes whose expression is associated with SCZ risk by TWAS is a powerful approach to prioritize SCZ risk genes. By integrating multiple published datasets from TWAS, gene coexpression, and differential gene expression analysis, we prioritized PCCB as a reliable SCZ risk gene. PCCB is a gene encoding β subunit of propionyl-coA carboxylase enzyme 28 , defect of which has been reported as a cause of propionic acidemia 33  To investigate PCCB's contribution to SCZ risk, we performed RNA-seq analysis and identi ed that PCCB knockdown in U2F hFOs affected expression of genes related to multiple neuronal functions and GABAergic synapse pathway. To con rm the RNA-seq results, we generated hFOs using another hiPSC line (ACS-1011) (Fig. S4A, B), nding that PCCB knockdown also led to decreased expression of GABA receptor genes, including GABRA1, GABRA2, GABRB2, and GABRB3 (Fig. S4C), as detected in U2F hFOs.
The downregulated GABAergic synapse pathway caused by PCCB knockdown attracted our attention, since GABAergic system dysfunction plays important roles in SCZ etiology 37,38 . We performed the metabolomic analysis and con rmed the decreased GABA levels in PCCB knockdown hFOs. The following electrophysiological analysis showed that PCCB knockdown led to hyper neuroactivity and decreased synchronization of neural network activities, cellular phenotypes reported to be associated with SCZ risk 32,39 . Through the hFOs-based multiomics analyses, we revealed the impacts of PCCB in neural functions and its connection to SCZ etiology, highlighting that PCCB may contribute to SCZ etiology through regulating the GABAergic system.
Since the GABA shunt connected the GABA metabolism pathway and TCA cycle 29,30 , we expected that PCCB regulates GABAergic system by affecting TCA cycle. As expected, PCCB knockdown lead to reduced production of SCOA and α-KG in TCA cycle. The α-KG produced from TCA cycle could be converted into SCOA or serve as a source for GABA synthesis 40 , the decrease of α-KG may be responsible for the reduced production of GABA. On the other hand, PCCB knockdown lead to the reduction of SCOA, which may exacerbate the entry of GABA into the TCA cycle through GABA shunt pathway and further reduced GABA content in the cytoplasm 40 . Overall, PCCB knockdown decreased GABA levels by reducing the content of α-KG and SCOA in TCA cycle (Fig. 4G). As mitochondrial dysfunction has been reported to be associated with GABA dysfunction 41 and etiology of SCZ 42,43 , our study provided evidence how mitochondrial dysfunction may contribute to SCZ risk.
In addition to mitochondrial dysfunction, one of the major effects of PCCB defect is the cellular accumulation of propanoic acid, propionyl carnitine, and other metabolites. Indeed, we did observe a dramatic increase of propanoic acid or propionyl carnitine in PCCB knockdown hFOs (Table S3). Interestingly, hFOs exposed to propanoic acid (3.5 uM) led to signi cantly decreased expression of GABA receptor genes (GABRA1, GABRA2, GABBR2, and GABBR3) (Fig. S5), which is consistent with those observed in PCCB knockdown hFOs. These results suggested that the accumulation of propanoic acid may mediate the effects of PCCB knockdown, and partially explain how propionic acidemia could lead to neuropathological symptoms. These results also suggested the potential effects of short-chain fatty acids on SCZ risk, since short-chain fatty acids including propanoic acid, acetic acid, and butyric acid were found to be upregulated in serum of SCZ patients 44 .
This study reveals the connection between PCCB and SCZ risk, some limitations also existed. First, PCCB knockdown affects multiple types of synapses, including GABAergic, glutamatergic, dopaminergic, and cholinergic synape, as revealed by RNA-seq analysis. The cell-type speci c effects of PCCB knockdown is unclear. Further investigations such as RNA-seq and other omics analysis at single-cell level are needed. Second, we showed that SCZ-associated SNP rs35874192, an eQTL SNP for PCCB, affected transcriptional activities. But whether SNP rs35874192 could affect PCCB expression in-vivo remains unclear. In the future, using CRISPR-Cas9 gene editing to con rm the regulatory effects of SNP rs35874192 on PCCB expression is needed.
In summary, this study used hFOs-based multi-omics analyses and revealed connection between PCCB and SCZ, highlighting that PCCB may contribute to SCZ etiology through regulating the GABAergic system and mitochondrial function.

Prioritization of SCZ risk genes and SNPs
We combined the published results from TWAS 8-11 , MR-JTI 15 , and SMR 2,10,12,13 analyses to prioritize SCZ risk genes with su cient supporting evidence. We also checked whether the prioritized genes are located in SCZ risk-associated gene coexpression modules or differentially expressed in postmortem brains of SCZ patients.
To prioritize SCZ risk SNPs, we rst collected top SNPs in the TWAS analysis. We then retrieved SNPs in LD (r 2 ≥ 0.6, European population genome) with the top SNPs. We prioritized candidate causal SNPs that likely affect gene expression in the brain using the following criteria: 1) candidate SNPs are eSNPs for the SCZ risk genes in the brain based on the BrainSeq 24 , GTEx 23 , or PsychENCODE eQTL data 45  and 1% penicillin/streptomycin (Gibco, 10378016).
As described in our previous study 46 , hNPCs were induced from U2F hiPSCs using the STEMdiff™ Neural Induction Medium (STEMCELL Technologies, 05835). The hNPCs were cultured in Matrigel-coated plate and maintained in the STEMdiff™ Neural Progenitor Medium (STEMCELL Technologies, 05833).
For DLRA, 1× 10 5 cells per well were plated into 24-well plate. After 24 hours of culture, 500 ng recombinant pGL3-basic luciferase reporter vector and 20 ng PRL-TK Renilla internal control vector for each well were co-transfected into cells using Lipofectamine™ 3000 (Invitrogen, L3000015). 36 hours post transfection, the Fire y and Renilla luciferase activities were measured on the LumiPro luminescence detector (Lu-2021-C001) using the DLRA kit (Vazyme, DL101-01). Experiments were conducted in three biological replicates.

Rt-qpcr Analysis
Total RNA was used to generate complementary DNA using HiScript III RT SuperMix for qPCR (+ gDNA wiper) (Vazyme, R323-01). RT-qPCR assay was performed using ChamQ SYBR qPCR Master Mix (Vazyme, Q711-02) on the Real-Time PCR System (Roche, LightCycler 480 II). GAPDH was used as the internal reference gene. At least three technical replicates were used in the RT-qPCR analysis. RT-qPCR primers are provided in Table S4.

Generation Of Hfos
The established PCCB knockdown and control hiPSCs were used to generate hFOs using the STEMdiff Dorsal Forebrain Organoid Differentiation Kit (STEMCELL Technologies, 08620) based on the manufacturer's instructions with some modi cations. Brie y, hiPSCs were dissociated into single cells using Accutase solution (Sigma-Aldrich, A6964) on day 0. 1x 10 4 cells per well were then plated into 96-well round-bottom ultra-low attachment plate (Corning, 7007) and fed with 50 µL Forebrain Organoid Formation Medium supplemented with 1x Penicillin/Streptomycin and 10 µM Y27632 (Selleck, SCM075). On day 3, each well was gently added with 50 µL fresh Forebrain Organoid Formation Medium without Y27632. On day 6, the medium was replaced with the Forebrain Organoid Expansion Medium. On day 25, the Forebrain Organoid Expansion Medium was replaced with the Forebrain Organoid Differentiation Medium. From day 43, organoids were cultured in the Forebrain Organoid Maintenance Medium. hFOs were characterized using immunostaining as described in our previous studies 46, 47 .

Bulk Organoid Rna-seq And Data Analysis
Two hFOs from PCCB knockdown or control group were randomly selected and pooled together as one mixed sample (day 60, N = 5 in each group) for total RNA extraction using the miRNeasy Mini Kit (Qiagen, 217004). RNA quality was evaluated on the Agilent 2100 Bioanalyzer system. RNA samples with RNA integrity numbers over 7 were used for RNA-seq (150 bp, paired-end) on the Illumina NovaSeq 6000 system.
Raw RNA-seq data were ltered to get clean reads using FastQC (v0.20.0). The clean reads were aligned to the human genome hg38 using STAR (v2.7.9a). Gene expression quanti cation was conducted using RSEM (v1.3.0) based on Gencode v40 comprehensive gene annotation. The lterByExpr function in the edgeR package (v3.36.0) was used to lter out low-expression genes. The sva function in the SVA package (v3.42.0) was used to estimate batch effect and other artifacts. Differential gene expression analysis between the PCCB knockdown and control group was performed using the DEseq2 package (v1.34.0). P values were adjusted using the Benjamini-Hochberg method. To annotate the functions of DEGs, we used the online tool WebGestalt2019 (http://www.webgestalt.org/) to perform GO and KEGG pathway enrichment analysis.

Ppi Analysis
The STRING database (v11.5) (http://www.string-db.org/) was used to construct a high-con dence (interaction score > 0.7) PPI network for the PCCB-induced DEGs. Active interaction sources included text-mining, experiments, databases, coexpression, neighborhood, gene fusion, and co-occurrence. The PPI network was visualized using Cytoscape (v 3.9.1). CytoHubba 49 , a Cytoscape plugin, was used to explore hub nodes in the PPI network.

Metabolomic Analysis Of Bulk Hfos
For metabolomic analysis, three hFOs (day 60) in PCCB knockdown or control group were randomly selected and pooled together as one mixed sample. Five mixed samples in each group were then used for HM400 metabolomic analysis (Beijing Genomics Institute, China). Brie y, hFOs or quality control samples were lysed in 140 µL 50% water/methanol solution. The lysate was centrifuged (18000g, 4℃, 20 min) to get the supernatant. The supernatant was used for derivatization reaction and then centrifuged at 4000g, 4℃, 10 min. The supernatant was further used for high performance liquid chromatography tandem mass spectrometer (LC-MS/MS) analysis on the SCIEX QTRAP 6500 + LC-MS/MS system. Parameters of liquid chromatographic column were BEH C18 (2.1 mm x 10 cm, 1.7um, Waters). The parameter of mass spectrometry was ESI+/ESI-. Content of metabolites (µmol/g) was quanti ed using the HMQuant software based on the formula (C*0.14/m), where C represents the calculated concentration (µmol/L) and m represents the sample weight (mg). Two-tailed t-test was used to identify differential metabolites between PCCB knockdown and control hFOs. Functional annotation for the differentially expressed metabolites was performed using the online tool MetaboAnalyst (https://www.metaboanalyst.ca/). A, B, and C DLRA results for PCCB eSNPs rs35874192 (A), rs900818 (B), and rs7349597 (C) in hNPCs and SH-SY5Y. The PRL-TK Renilla vector was used as internal control. Data are shown as Mean ± SEM. Unpaired two tailed t test was used for comparison between two groups. **P < 0.01, ***P < 0.001. Reference eQTL plots in this gure were downloaded from the GTEx portal (https://www.gtexportal.org/home/) and BrainSeq phaseIeQTL data (http://eqtl.brainseq.org/phase1/eqtl/).

Figure 2
Functional effects of PCCB knockdown in hFOs.
A Work ow of hFOs culture in this study. Scale bar, 200 μm and 300 μm. B, C RT-qPCR analysis for PCCB expression in hiPSCs (B) and hFOs (C). D Immunostaining characterization for hFOs. Cortical plate marker, MAP2; intermediate zone marker, TBR2; ventricular zone markers, SOX2 and ki67; forebrain-speci c markers, FOXG1 and PAX6. Scale bar, 50 μm. EPCA plot for the PCCB knockdown and control hFOs. PCCB-NC was shown with green dots and PCCB-G1 was shown with red dots. F Volcano plot of DEGs between the PCCB knockdown and control hFOs. Upregulated genes are shown with red dots and downregulated genes are shown with blue dots. G, H GO and KEGG analysis for PCCBinduced upregulated (G) and downregulated DEGs (H) respectively. GO terms are shown with red bars, and KEGG are shown with blue bars. I PPI network analysis for 350 shared PCCB-induced downregulated genes between PCCB-G1 and PCCB-G2 hFOs. The top 10 hub nodes are shown with orange nodes. J RT-qPCR analysis for GABRA1, GABRA2, GABRB2, and GABRB3 (The hub nodes in the PPI network). The twotailed student's t-test was used to assess difference between the PCCB-NC and PCCB-G1 or PCCB-G2 group. **P< 0.01, ***P < 0.001, ****P < 0.0001. PCCB knockdown leads to mitochondrial dysfunction in hFOs.
A Bright elds of hFOs cultured in 24-well MEA plate. Scale bar, 100 μm. B Representative burst traces for individual electrodes recorded from hFOs. C Schematic diagram of a single unit event. D Raster plot diagram of synchronized burst activity. Each pink box represents a synchronized burst. E, F, and G PCCB knockdown increased the number of spikes (E) and mean neuron ring rate (F), but reduced synchronized burst activity in hFOs (G). Data were shown with Mean ± SEM (averaged organoids N = 8). The two-tailed student's t-test was used to assess difference between the PCCB-NC and PCCB-G1 or PCCB-G2 group. *P< 0.05, ****P < 0.0001.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.