Activity-Dependent Changes in Gene Expression in Schizophrenia Human-Induced Pluripotent Stem Cell Neurons

Key Points

Question  Is the neuronal activity-dependent change of gene expression dysregulated in schizophrenia?

Findings  In this genetic analysis, candidate genes converged to gene networks that are associated with differential effect of activity–dependent changes of gene expression in schizophrenia.

Meaning  Candidate genes for schizophrenia might act through dysregulation of the ability of neurons to activate molecular processes in response to depolarization.

Importance  Schizophrenia candidate genes participate in common molecular pathways that are regulated by activity–dependent changes in neurons. One important next step is to further our understanding on the role of activity-dependent changes of gene expression in the etiopathogenesis of schizophrenia.

Objective  To examine whether neuronal activity-dependent changes of gene expression are dysregulated in schizophrenia.

Design, Setting, and Participants  Neurons differentiated from human-induced pluripotent stem cells derived from 4 individuals with schizophrenia and 4 unaffected control individuals were depolarized using potassium chloride. RNA was extracted followed by genome-wide profiling of the transcriptome. Neurons were planted on June 21, 2013, and harvested on August 2, 2013.

Main Outcomes and Measures  We performed differential expression analysis and gene coexpression analysis to identify activity-dependent or disease-specific changes of the transcriptome. Gene expression differences were assessed with linear models. Furthermore, we used gene set analyses to identify coexpressed modules that are enriched for schizophrenia risk genes.

Results  We identified 1669 genes that were significantly different in schizophrenia-associated vs control human-induced pluripotent stem cell–derived neurons and 1199 genes that are altered in these cells in response to depolarization (linear models at false discovery rate ≤0.05). The effect of activity-dependent changes of gene expression in schizophrenia-associated neurons (59 significant genes at false discovery rate ≤0.05) was attenuated compared with control samples (594 significant genes at false discovery rate ≤0.05). Using gene coexpression analysis, we identified 2 modules (turquoise and brown) that were associated with diagnosis status and 2 modules (yellow and green) that were associated with depolarization at a false discovery rate of ≤0.05. For 3 of the 4 modules, we found enrichment with schizophrenia-associated variants: brown (χ² = 20.68; P = .002), turquoise (χ² = 12.95; P = .04), and yellow (χ² = 15.34; P = .02).

Conclusions and Relevance  In this analysis, candidate genes clustered within gene networks that were associated with a blunted effect of activity-dependent changes of gene expression in schizophrenia-associated neurons. Overall, these findings link schizophrenia candidate genes with specific molecular functions in neurons, which could be used to examine underlying mechanisms and therapeutic interventions related to schizophrenia.

Genetic studies have implicated numerous risk variants for schizophrenia.^1-7 One of the next challenges is to further understand the biological mechanisms of the large number and diversity of genes that are associated with schizophrenia. To do that, we need to generate functional data capturing molecular pathways that are relevant to schizophrenia and examine the involvement of multiple candidate genes. Such strategy holds the potential to advance our understanding of the underlying molecular basis of schizophrenia and a way to develop novel treatments.

Studies suggest that a convergent molecular pathway dysregulated in schizophrenia is the signaling network that modulates synaptic strength or the network of genes that are targets of fragile X mental retardation protein.^2,7 In addition, these networks are enriched for genes affected by mutations in autism and intellectual disability,⁷ providing further insight for a pathophysiology shared among neurodevelopmental disorders. Interestingly, both gene networks are composed of many proteins for which expression is regulated by neuronal activity,⁸ further suggesting that schizophrenia might be driven by dysregulation of neuronal activity–dependent synapse development and function.

Neurons differentiated from human-induced pluripotent stem cells (hiPSCs) can assess human brain cellular properties⁹ and have the potential to detect changes of gene expression in response to neuronal depolarization. We used hiPSC neurons derived from patients with schizophrenia and unaffected control individuals and identified disease-specific differences in gene expression as well as blunted effects of activity-dependent changes of gene expression in schizophrenia compared with control samples. Furthermore, these differentially expressed genes are coexpressed in modules that are highly enriched for genes affected by genetic risk variants in schizophrenia and other neurodevelopmental disorders. In this study, we aimed to determine whether schizophrenia candidate genes converge to gene networks that are associated with the differential effect of activity-dependent changes of gene expression in schizophrenia.

Forebrain neural progenitor cells (NPCs) previously differentiated from case and control hiPSCs as reported^9,10 were subjected to 6 weeks of neuronal differentiation and maturation (eAppendix in the Supplement). Because all research described herein was performed on deidentified human samples obtained for broadly consented scientific research by either American Type Culture Collection or the Coriell Cell Repository, it was found to be exempt by the internal review committee of the Icahn School of Medicine at Mount Sinai. Neurons were treated with 50mM potassium chloride (or phosphate-buffered sulfate vehicle control) for the final 3 hours prior to harvest. Five hundred nanograms of total RNA were used as input material for library preparation using the TruSeq Stranded Total RNA Kit (Illumina).

Reads were mapped to hg19 reference genome using TopHat. After exploratory analysis, we identified and included the following 5 covariates in the differential gene expression analysis: diagnosis, treatment, sex, age, and RNA integrity number. For each transcript, we fit linear regression models for the effect on gene expression of each variable, using the limma package: gene expression ~ diagnosis + treatment + covariates.

Differentially expressed genes were determined at false discovery rate (FDR)¹¹ of 0.05 or less. The χ² results were determined using the Fisher method.

We validated transcripts that were significant for disease status or altered after depolarization using real-time polymerase chain reaction. Briefly, 500 ng of total RNA were used for complementary DNA (cDNA) synthesis, and each cDNA sample was amplified in triplicate using SYBR Green (ThermoFisher Scientific). The primers used for this analysis are described in eTable 1 in the Supplement.

We constructed unsigned gene coexpression networks using the weighted gene coexpression analysis (WGCNA) package in R (R Programming), starting with the normalized and residualized (removing effect of age, sex, and RNA integrity number) expression data for 13 903 genes. For these data, we used an R² cutoff of 0.85, which corresponded to a selection of β = 7. Ordered from largest (the module containing the most genes) to smallest, each module was sequentially assigned a color name. The less well-connected genes were arbitrarily grouped in the gray module.

Enrichment analysis was performed using gene sets for known molecular pathways and biological processes, including gene ontology (http://www.geneontology.org)¹²; the Reactome data set (http://www.reactome.org)¹³; and the HUGO Gene Nomenclature Committee gene families (http://www.genenames.org).¹⁴ In addition, to further characterize the data-driven coexpressed modules, we generated a group of gene sets derived from previous schizophrenia genetic findings, cell type–specific studies, or coexpression analyses. The genes in each module were tested for overlap using Fisher exact test and FDR correction across all modules and all gene sets tested.

We tested the enrichment of modules for genes found in genetic loci previously associated with schizophrenia, including (1) 108 loci discovered in a common variant genome-wide association study¹; (2) a literature consensus of 12 copy number variant regions collated from numerous rare copy number variant studies¹⁵; and (3) 756 nonsynonymous (mostly missense but also including 114 loss-of-function [nonsense, essential splice site, or frameshifting indels]) de novo mutations.^2-6 In addition, we used published de novo mutations across 3 neurodevelopmental disorders, including autism spectrum disorders,^16-19 intellectual disability,^20-23 and epilepsy.^24,25

In previous publications, we directly reprogrammed fibroblasts from 4 patients with schizophrenia and 6 control individuals into hiPSCs and differentiated these disorder-specific hiPSCs into forebrain NPCs¹⁰ and neurons.⁹ As previously reported, spatial patterning, replication, and propensity toward neuronal differentiation of NPCs does not appear to vary based on psychiatric diagnosis.^9,10 When differentiated to neurons, forebrain NPCs yield a population of forebrain hiPSC neurons that are vesicular glutamate transporter 1–positive and so are presumably excitatory glutamatergic neurons, although approximately 30% of neurons are glutamate decarboxylase 67–positive (GABAergic).⁹

In this study, we compared global transcription of hiPSC forebrain neurons from 4 control participants and 4 patients with schizophrenia (eTable 2 in the Supplement) with or without potassium chloride-induced depolarization (PCID). RNA sequencing was performed and we obtained, on average, 21.5 million paired-end reads, among which 94% were mapped on the genome (eTable 2 in the Supplement). Following data normalization, there were 13 903 genes for analysis, of which 12 678 were protein coding. Multidimensional scaling separated control samples from schizophrenia samples and phosphate-buffered sulfate from potassium chloride–treated samples (eFigure 1 in the Supplement). We first explored transcriptional changes related to schizophrenia by identifying differentially expressed genes (DEGs) between schizophrenia and control neuronal samples. We identified 1669 DEGs at FDR of 0.05 or less (eTable 3 in the Supplement), as illustrated by the heat map (Figure 1A) and volcano plots (Figure 1B). Among the 1669 DEGs, 854 were upregulated and 815 were downregulated in schizophrenia with a moderate to strong effect (mean log₂ fold change [FC], 3.91; range, 0.77-14.50). eFigure 2 in the Supplement shows illustrative examples for 2 DEGs that have also been previously associated with schizophrenia genetic loci¹: CACN2B (log₂ FC, −2.1; P = 8.7 × 10⁻⁴ at FDR 1.6%) and GRM3 (log₂ FC, 3.3; P = 2.5 × 10⁻³ at FDR 3.0%). The results from the DEG analysis were validated for 4 transcripts using quantitative polymerase chain reaction (eFigure 3 in the Supplement).

We examined differential expression in a previous sample that generated gene expression data using microarrays from 4 patients with schizophrenia (12 microarrays) and 4 control individuals (12 microarrays).⁹ While these arrays differed from RNAseq in their capture features, the correlation of test statistics for differential expression in this current data set compared with the previous was significant: Spearman correlation ρ = 0.62 (P = 5.7 × 10⁻¹⁶¹) for the subset of significantly DEGs also present in the previous data set (n = 1509 genes) (eFigure 4 in the Supplement). Pathway enrichment analysis across a broad set of pathways was conducted for interpretation of the list of DEGs. The most enriched categories included pathways related to organization of extracellular matrix, cell adhesion, and binding of calcium ion (eFigure 5 in the Supplement). We note that pathway enrichment analysis using the DEGs from the microarray analysis show enrichment for similar pathways, with the exception of the biological process related to the response to retinoic acid that shows higher enrichment in this data set (data not shown). Overall, these results indicate that only a few specific pathways are affected in schizophrenia vs control neurons.

We next examined transcriptional changes related to PCID. We identified 1199 DEGs at FDR of 0.05 or less (eTable 4 in the Supplement). Among the 1199 genes, 595 were upregulated and 604 were downregulated in response to PCID with a moderate effect (average log₂ FC, 1.26; range, 0.45-4.82) (Figure 2A). Unsupervised clustering shows a clear separation of PCID signatures in controls but less clear in schizophrenia, indicating a less significant effect of PCID on gene expression in schizophrenia (Figure 2B). eFigure 6 in the Supplement shows 2 illustrative examples for SLC38A2 (log₂ FC, 4.2; P = 1.4 × 10⁻⁸ at FDR 0.007%) and THNSL1 (log2 FC, 2.7; P = 2.1 × 10⁻⁷ at FDR 0.02%) that are upregulated and downregulated in response to PCID, respectively.

We examined PCID signatures in a previous sample that generated gene expression data in mouse cortical neurons.²⁶ We found a moderate and significant correlation between our PCID signatures and mouse cortical neuron signatures after PCID for 1 hour (Spearman correlation ρ = 0.23; P = 7.7 × 10⁻¹⁵) or 6 hours (Spearman correlation ρ = 0.27; P = 2.6 × 10⁻²⁰) (eFigure 7 in the Supplement). Pathway enrichment analysis showed enrichment for molecular functions and biological processes related to transcription factor activity, regulation of transcription, and signaling by ERBB4 and NRG (eFigure 8 in the Supplement).

Based on the more prominent separation of PCID signatures in control samples compared with schizophrenia, we performed separate DEG analysis in each group (Figure 2B). Remarkably, whereas 594 genes were PCID DEGs in controls at FDR of 0.05 or less (eTable 5 in the Supplement), only 59 genes were differentially expressed in schizophrenia (eTable 6 in the Supplement). This was not simply an issue of statistical thresholds because relaxing the statistical criteria for differential expression (FDR ≤0.1) identified 1112 differentially expressed genes in control participants and only 169 in patients with schizophrenia, confirming the large difference observed in PCID signatures. We note that PCID signatures were consistent in terms of directionality (upregulated or downregulated) in both control participants and patients with schizophrenia (Figure 3A). The results from the DEG analysis were validated for 4 transcripts using quantitative polymerase chain reaction (eFigure 9 in the Supplement).

To further quantify the disease effect (control and schizophrenia) on the PCID transcriptional variability, we applied a linear mixed model for each gene and quantified the total variance of PCID within each diagnostic group.²⁷ Genome-wide variation across PCID accounted for a median of 7.1% variation in control samples. In schizophrenia samples, the variation explained was 4.9%, confirming our previous observation of attenuated effect of PCID on gene expression in schizophrenia compared with controls (Figure 3B). We then ranked genes based on the difference of the variance of PCID signature in control participants vs patients with schizophrenia. We identified 2040 genes (control signatures) where more than 20% and less than 20% of their variation was explained by treatment in control participants and patients with schizophrenia, respectively. On the other hand, only 1111 genes (schizophrenia signatures) had less than 20% and more than 20% of their variation explained by treatment in control participants and patients with schizophrenia, respectively. The top 20 control signatures and an illustrative example for KIF5B gene are provided (eFigures 10 and 11 in the Supplement). Pathway enrichment analysis on the top 5% of control signatures showed enrichment for molecular functions and biological processes related to transcription factor activity and regulation of transcription (eFigure 12 in the Supplement).

To identify discrete groups of coexpressed genes showing transcriptional differences relevant to diagnosis or PCID, we constructed a coexpression network using the entire data set. Genes were clustered into 15 coexpression modules based on high similarity of expression patterns across samples. Genes that clustered into specific modules based on similar coexpression patterns were also DEGs for diagnosis or PCID status, indicating that they participate in common biological processes (Figure 4A).

We then examined whether the identified coexpression networks recapitulate molecular processes related to diagnosis and PCID status. We used the module eigengene (the first principal component of the expression pattern of the corresponding module) to summarize gene expression trajectories across samples and evaluated the relationship of the 15 module eigengenes with diagnosis and PCID status (Figure 4B; eTable 7 in the Supplement). We found 2 modules strongly associated with diagnosis status at FDR of 0.05 or less. One module was upregulated in schizophrenia (turquoise) and was enriched for ribosome markers and for genes belonging to 3′ untranslated region–mediated translational regulation (eFigure 13 in the Supplement). The other module was downregulated in schizophrenia (brown) and was enriched for neuronal markers, as well as for genes belonging to fragile X mental retardation protein targets (eFigure 13 in the Supplement).

We also found 2 modules strongly associated with PCID status at FDR of 0.05 or less. One module was upregulated with potassium chloride treatment (yellow) and was enriched for genes belonging to tRNA processing and ZNF gene family (eFigure 13 in the Supplement). The other module was downregulated with potassium chloride treatment (green) and was enriched for genes belonging to regulation of transcription and ZNF gene family (eFigure 13 in the Supplement). Both modules were significant when we examined the association of module eigengene with PCID in control samples, while only the green module was significant in schizophrenia. These results are consistent with the less significant effect of PCID on gene expression (and module eigengene) in schizophrenia samples compared with control samples.

To further interpret the coexpressed modules, we examined whether modules are enriched for genes that have been previously associated with schizophrenia and other neurodevelopmental illnesses (autism, epilepsy, and intellectual disability) (Figure 5; eTable 8 in the Supplement). Common schizophrenia risk variants defined based on the most recent large-scale genome-wide association study from Psychiatric Genomics Consortium 2¹ were enriched for 5 modules at P < .05. Interestingly, 3 of the 5 modules were associated with disease status: brown (observed to expected ratio, 1.51; P = .003) and turquoise (observed to expected ratio, 1.30; P = .02) or PCID-yellow (P = .03). A significant overlap of the yellow module (observed to expected ratio, 1.51; P = .05) was identified with genes that lie within genomic regions that have been associated with large cytogenomic deletions in schizophrenia; no significant effect was observed for schizophrenia-associated duplicated regions.

The brown module was significantly enriched for nonsynonymous and loss-of-function de novo mutations across multiple neurodevelopmental diseases, including autism, epilepsy, and intellectual disability, but not schizophrenia. In addition, the turquoise module was significantly enriched for de novo mutations in autism. Schizophrenia de novo mutations were enriched for 3 modules; none of these modules were associated with disease or PCID status. Finally, we examined the effect of schizophrenia-associated genetic signal in aggregate and identified a significant effect for 3 of the 4 modules that were associated with disease or PCID status: brown (χ² = 20.68; P = .002), turquoise (χ² = 12.95; P = .044), and yellow (χ² = 15.34; P = .0218). The hub nodes (most connected genes) are illustrated for the brown, turquoise, and yellow modules (eFigure 14 in the Supplement).

Neuronal depolarization regulates synaptic function by inducing gene expression and protein synthesis of synaptic molecules in dendrites.⁸ A previous study reported activity-dependent changes in neurotransmitter release in neurons differentiated from these same hiPSC lines.²⁸ Here, we report changes in gene expression in response to depolarization. These genes participate in molecular functions and biological processes related to regulation of transcription, and they show consistent changes with gene expression studies in mouse cortical neurons.²⁶ Thus, the transcriptome of hiPSC neurons is induced in response to neuronal depolarization, indicating that they are a relevant platform for the study of activity-dependent changes in neuropsychiatric diseases.

Compared with control samples, hiPSC neurons derived from patients with schizophrenia showed alterations of gene expression consistent with a previous study⁹ and are enriched in molecular pathways related to organization of extracellular matrix and cell adhesion. These extracellular matrix proteins include gene families, such as neural cell adhesion molecules, neuroligins, and neurexins, that are critical for cellular differentiation and migration²⁹ and are postulated as a possible mechanism for the etiopathogenesis of schizophrenia.³⁰ Previous studies have demonstrated that schizophrenia NPCs have abnormal expression for extracellular matrix and cell adhesion genes,^10,31 accompanied by deficits in migration.¹⁰

By examining the higher-order organization of the transcriptome, we identified specific modules to be associated with disease status and PCID. We identified 2 modules that were significantly deregulated in schizophrenia and 2 modules that were altered in response to PCID. Similarly, to differential expression analysis, we observed a blunted effect of activity-dependent changes of module activity in schizophrenia. By testing the enrichment of modules with schizophrenia candidate genes, we found a strong enrichment of genetic variants with the disease-associated modules and the PCID-associated network that was only activated in control samples but not schizophrenia. These data support the notion that the impairment of different but not mutually exclusive molecular pathways contributes to the etiology of schizophrenia. For instance, schizophrenia may arise from early defects in neurodevelopment, followed by global dysregulation of synaptic function and impairment of neurotransmission. This is consistent with findings in the brown module, which participates in neurodevelopment and synaptic function, has a strong enrichment for schizophrenia variants, and is associated with the disease status. Alternatively, another possibility is that causative factors involve dysregulation of activity-dependent signaling pathways either locally at the synapse or globally, leading to deficits in activity-dependent gene transcription. This is consistent with results in the yellow module, which is enriched for schizophrenia risk variants and has a strong association with PCID in controls.

While our study shows convergence of genetic variants to specific dysregulated coexpressed networks, it fails to identify driver genes. This can be the focus of future studies, where incorporating larger sample sizes and including additional times for PCID treatment would allow us to generate directed networks, such as causal probabilistic networks, and identify specific genes that drive abnormalities in schizophrenia. Finally, future studies should collect additional functional data, such as quantification of secreted neurotransmitters and electrophysiological measurements, which would allow investigating the effect of dysregulated gene expression on neurons’ physiology.

In this study, we examined the effect of depolarization on gene expression in hiPSC neurons from patients with schizophrenia and control individuals. Schizophrenia neurons failed to induce the gene expression of specific transcripts to the same degree that occured in control samples. These transcripts cluster in coexpressed networks that contain candidate genes for schizophrenia and other neurodevelopmental diseases. Overall, these findings provide further insight into the function and mechanism of schizophrenia candidate genes.

Corresponding Author: Panos Roussos, MD, PhD, Departments of Psychiatry and Genetics and Genomic Science, Icahn School of Medicine at Mount Sinai, 1470 Madison Ave, New York, NY 10029 (panagiotis.roussos@mssm.edu).

Accepted for Publication: August 21, 2016.

Published Online: October 12, 2016. doi:10.1001/jamapsychiatry.2016.2575

Author Contributions: Drs Roussos and Brennand had full access to all of the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis. Drs Barry and Brennand contributed equally.

Concept and design: Roussos, Barry, Brennand.

Acquisition, analysis, or interpretation of data: All authors.

Drafting of the manuscript: Roussos.

Critical revision of the manuscript for important intellectual content: All authors.

Statistical analysis: Roussos, Guennewig.

Administrative, technical, or material support: Roussos, Kaczorowski, Brennand.

Study supervision: Roussos.

Conflict of Interest Disclosures: None reported.

Funding/Support: This work was supported by the National Institutes of Health grants R01AG050986 and R01MH109677 (Dr Roussos) and R01MH101454 and R01MH106056 (Dr Brennand), New York Stem Cell Foundation (Dr Brennand), Brain Behavior Research Foundation (Drs Roussos and Brennand), Alzheimer's Association grant NIRG-340998 (Dr Roussos), the Veterans Affairs merit grant BX002395 (Dr Roussos), the Swiss National Science Foundation grant P2EZP3_152143 (Dr Guennewig), and a Geoff and Dawn Dixon fellowship (Dr Guennewig). This work was supported in part through the computational resources and staff expertise provided by Scientific Computing at the Icahn School of Medicine at Mount Sinai. Research reported in this article was supported by the Office of Research Infrastructure of the National Institutes of Health award S10OD018522.

Role of the Funder/Sponsor: The funding sources had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Disclaimer: The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.