Co-expression network of neural-differentiation genes shows specific pattern in schizophrenia

Background Schizophrenia is a neurodevelopmental disorder with genetic and environmental factors contributing to its pathogenesis, although the mechanism is unknown due to the difficulties in accessing diseased tissue during human neurodevelopment. The aim of this study was to find neuronal differentiation genes disrupted in schizophrenia and to evaluate those genes in post-mortem brain tissues from schizophrenia cases and controls. Methods We analyzed differentially expressed genes (DEG), copy number variation (CNV) and differential methylation in human induced pluripotent stem cells (hiPSC) derived from fibroblasts from one control and one schizophrenia patient and further differentiated into neuron (NPC). Expression of the DEG were analyzed with microarrays of post-mortem brain tissue (frontal cortex) cohort of 29 schizophrenia cases and 30 controls. A Weighted Gene Co-expression Network Analysis (WGCNA) using the DEG was used to detect clusters of co-expressed genes that werenon-conserved between adult cases and controls brain samples. Results We identified methylation alterations potentially involved with neuronal differentiation in schizophrenia, which displayed an over-representation of genes related to chromatin remodeling complex (adjP = 0.04). We found 228 DEG associated with neuronal differentiation. These genes were involved with metabolic processes, signal transduction, nervous system development, regulation of neurogenesis and neuronal differentiation. Between adult brain samples from cases and controls there were 233 DEG, with only four genes overlapping with the 228 DEG, probably because we compared single cell to tissue bulks and more importantly, the cells were at different stages of development. The comparison of the co-expressed network of the 228 genes in adult brain samples between cases and controls revealed a less conserved module enriched for genes associated with oxidative stress and negative regulation of cell differentiation. Conclusion This study supports the relevance of using cellular approaches to dissect molecular aspects of neurogenesis with impact in the schizophrenic brain. We showed that, although generated by different approaches, both sets of DEG associated to schizophrenia were involved with neocortical development. The results add to the hypothesis that critical metabolic changes may be occurring during early neurodevelopment influencing faulty development of the brain and potentially contributing to further vulnerability to the illness. Electronic supplementary material The online version of this article (doi:10.1186/s12920-015-0098-9) contains supplementary material, which is available to authorized users.


Background
Schizophrenia is a severe and debilitating neuropsychiatric disorder with a worldwide prevalence of approximately 1 % and an estimated heritability of 80-85 % [1]. Multiple factors contribute to its pathogenesis, including genetic and environmental factors [2][3][4], which are involved in the disturbance of the neuronal network formation and brain function [5].
Although the genetic component involvement is clear, the precise molecular and cellular defects that precedes illness onset and the sequence of mechanisms that initiate or attenuate disease progression, particularly in neurons, are still unknown. In addition to findings from animal models [6][7][8], post-mortem human brain studies in schizophrenia revealed significant alterations in genes with high expression in multiple cortical and hippocampal regions [9][10][11][12][13]. Data from prefrontal cortex obtained post-mortem revealed 28 modules of co-expressed genes associated to oxidative phosphorylation, energy production and metabolism, post-translation modifications, neurogenesis and neuronal differentiation. Additionally, these modules showed a robust effect of age suggesting that genes related to nervous system development process and neuronal differentiation are deregulated in schizophrenia [14]. A co-expression module analysis from patients with schizophrenia (SCZP) and controls of seven studies showed that the top-modules of each study displayed a significant enrichment of neuronal markers as well as genes related to the electron transport chain and oxidative phosphorylation biological processes [15]. Eventually, these alterations, through the analyses of genes and proteins, can be found in the peripheral blood suggesting that they are potential molecular markers for diagnosis or follow-up [16][17][18][19].
Promising approaches involve the study of tissue from the olfactory system by nose biopsies [20] and the development of human induced pluripotent cells (hiPSC) techniques. The later demonstrated that functional primarily glutamatergic neurons differentiated from hiPSC were able to fire action potentials after eight weeks in culture, showing the possibility of studying schizophrenia through this technique [21]. SCZP-hiPSC-derived neurons display diminished neuronal connectivity with reduced neuritis and changes in Wnt and cAMP signaling pathways when compared to controls [22]. Additionally, RNA sequencing analysis of hiPSCs-derived neurons highlighted the importance of long non-coding RNAs in neurogenesis [23].
In a previous study, we showed that neural progenitor cells (NPCs), derived from hiPSCs from skin fibroblasts of a schizophrenic patient had a 2-fold increase in extramitochondrial oxygen consumption compared to normal controls, and correspondingly elevated levels of reactive oxygen species (ROS), which was reduced to control levels by addition of valproic acid [24]. Mitochondrial dysfunction was also associated with impaired differentiation of neurons into both mature dopaminergic and glutamatergic phenotypes in schizophrenia-derived hiPSCs [25].
Although genes involved with ROS and other molecular mechanisms in hiPSC and NPC have been studied in patients with schizophrenia [26], the ones specifically disrupted in neuronal differentiation were not assessed directly in adult brain samples. The aim of this study was to disclose gene expression changes during neuronal differentiation disrupted in schizophrenia and to expand these findings to the brain tissue in adults. To achieve that, genes identified as involved with neuronal differentiation exclusively in schizophrenia were analyzed by coexpression network in post-mortem brain samples derived from adult patients and controls. Then, we identified the biological processes associated with a less conserved module between cases and controls brain samples.

Tissue collection
In short, biopsies of primary human fibroblasts from a 48 year-old woman with DSM-IV-TR schizophrenia defined as treatment-resistant and under clozapine treatment (SZCP) were collected in parallel with a control subject (CON) that negative for any major lifetime DSM-IV-TR diagnosis. Fibroblasts underwent primary culture procedures to generate hiPSCs and were subsequently differentiated into neurons (NPC). Details were previously described [24].
Total RNA was isolated from three plates of cell culture (with molecular confirmation for the cell type that was being studied [24]) using Trizol (TRI Reagent, Sigma-Aldrich, St Louis, MO, USA) according to the manufacturer's instructions. RNA quantity and integrity were assessed by spectrophotometry (NanoDrop ND-2000 UV-Vis Spectrophotometer, NanoDrop Technologies) and microfluidics-based electrophoresis (Agilent 2100 Bioanalyzer; Agilent Technologies, Santa Clara, California), respectively, and obtained OD of~2.0 and RIN > 9. Genomic DNA was extracted using standard proteinase K digestion followed by phenol/chloroform extraction. DNA quantity and integrity were assessed by spectrophotometry (NanoDrop ND-2000 UV Spectrophotometer, NanoDrop Technologies) and agarosis gel 0.8 %.
The study protocol was approved by the ethics committee of the Hospital de Clínicas de Porto Alegre, Porto Alegre (RS) with all subjects providing written informed consent co-signed by close relatives (Register 09-108). The study was performed in accordance with the Declaration of Helsinki.

Comparative genomic hybridization based on microarray (array-CGH)
We performed comparative genomic hybridization based on microarrays (array-CGH) in a commercial whole-genome 180 K platform containing 180,000 oligonucleotide probes (design 30864, Agilent Technologies, Palo Alto, California) using as reference DNA a commercially available human pool of healthy individuals (Promega). Purification, cohybridization of labelled test and reference DNA samples and washing were carried out according to the manufacturer's instructions. Scanned images were processed using the Feature Extraction 10.7.3.1 software (Agilent Technologies, Palo Alto, California).

Methylation experiments
Bisulfite-conversion of 500 ng of DNA was performed according to the manufacturer's recommendations (EZ DNA Methylation Gold Kit, Zymo Research, Irvine, California) for the Illumina Infinium Assay. The incubation profile was 16 cycles at 95°C for 30 s, 50°C for 60 min and a final holding step at 4°C for 10 min. In total, 10 μl of bisulfite-converted DNA were used for hybridization on Infinium HumanMethylation 450 BeadChip with 485,557 probes, following the Illumina Infinium HD Methylation protocol (Illumina, Inc, San Diego, California). All arrays passed Illumina standard quality control metrics and were included in this study.

Microarray experiments
For mRNA expression, samples were prepared according to the Agilent mRNA Microarray System protocol using 250 ng of total RNA for amplification and labelled using the Agilent Low RNA Input Fluorescent Linear Amplification Kit. Samples were hybridized in the Agilent Whole Human Genome Microarray Kit, 4 × 44 K-G4112F in the Agilent SureHyb chambers (Agilent Technologies, Palo Alto, California) for 20 h at 55°C, and the platforms were washed with the manufacturer's washing buffers. Controls were used in all labelling reactions. Scanned images of the arrays were processed using Feature Extraction 10.7.3.1 software (Agilent Technologies, Palo Alto, California) with default parameters. Quality control of the probes (i.e. low intensity, saturation, controls, etc.) selected 30,635 probes corresponding to 22,186 transcripts for the analyses. From those, 14,491 probes correspond to unique genes and 7695 were not associates with an ENTREZID. Four genes had more than one ENTREZID with the same Official Symbol and were excluded from the analyses.
Frontal cortical samples from a previous study were used for further analyses, all presenting good quality total and amplified RNA [27]. Postmortem brain tissue was donated by The Stanley Medical Research Institute's brain collection. Details about ethics regarding these samples are available at the Neuropathology Consortium (http://www.stanleyresearch.org). Comparing patients with schizophrenia and controls, there were no significant differences concerning age, brain weight or PMI, with only pH showing a significant difference (p = 0.01). In short, 30 controls subjects (CTS) and 29 patients with schizophrenia (SZP) from the Stanley Foundation were hybridized in the Agilent 4x44K human oligonucleotide microarray (Agilent 4112 F, Agilent Technologies, Palo Alto, California) following manufacturer procedures. For analysis of replicate spots, the average intensity after background correction was calculated, followed by normalization using locally weighted linear regression (LOWESS) with α = 0.2 within slides using the R software version 2.11.1 (R Development Core Team, 2010).
Data were deposited in NCBI Gene Expression Omnibus and are accessible through GEO Series accession number GSE62191 and GSE62105.

Data analyses
DNA copy number variations (CNVs) were identified with Nexus software 7.0 (Biodiscovery, Hawthorne, CA, USA) using the FASST2 segmentation algorithm. To be considered a CNV, the criteria were: a minimum of 3 consecutive affected probes (effective resolution of 50 Kb for CNV calling), a significance threshold set at 5.0e −8 , and threshold log 2 ratio Cy3/Cy5 of 0.4 for gain and −0.4 for loss. Data from sex chromosomes were excluded. Gene annotation was performed according to the GRCh37/hg19using the University of California Santa Cruz Genome Browser (http://genome.ucsc.edu/) and the Database of Genomic Variants (DGV; May 2014 update).
Due to the low number of samples, we simplified the methylation analyses. CpGs with β-values ≤0.2 were considered hypo-methylated and with β-values ≥0.8 were considered hyper-methylated, which were used to construct Venn-diagrams to compare the data in fibroblasts, hiPSC and NPC.
For gene expression analyses, it was applied HT-self [28]. Self-self experiments were performed by labelling the hIPSC from SZCP with either Cy3 or Cy5 dyes and hybridizing simultaneously on the same microarray slide to determine the intensity dependent cut-offs. These cut-offs were applied to the non-self-self experiment that was performed in triplicate (SCZP × CON). To identify differentially expressed genes, three criteria were applied: a. each gene must be represented by at least two probes; b. more than 50 % of the probes representing that gene presented reliable signal; c. there was at least 80 % of signal agreement among probes. To avoid technical variability, only genes that reached the three criteria in all replicates were selected.
To investigate if genes identified in this study could have a role during foetal development of the prefrontal cortex (dorsolateral prefrontal cortex, ventrolateral prefrontal cortex, anterior (rostral) cingulate (medial prefrontal) cortex), they were searched in the BrainSpan database (Allan Brain Atlas) [29]. Genes with at least 1 RPKM in more than 75 % of the samples or those who had an expression interquartile range >0.5 were considered for the analyses. Genes were assigned as foetal or non-foetal according to the development stage of the sample. It was considered differentially expressed those genes with fold-change >2.
The Wilcoxon Rank Sum Test was applied using TMEV software package [30] to identify differentially expressed genes between 29 SZP and 30 CTS considering those with p < 0.01 as significantly differentially expressed. As pH was significantly different in the brains between cases and controls, we accessed the DEG which expression was possibly being affected by pH. Hence, a multiple regression analysis was applied for all expressed genes in the microarray platform with a further comparison to 1000 permutated datasets. We observed that the difference of pH in the brain affects around 20 % of expressed genes (F, p < 0.05) of the whole array. These genes were compared to the DEGs between SZP and CTS.
The DEG sets found by both analyses were compared to the findings of a published study that defined 17 modules comprised of co-expressed genes associated with human cortical development. Briefly, to define human brain development modules, RNA-sequencing data from neocortical regions of the BrainSpan were analyzed through co-expression networks using the R package Weighted Gene Co-expression Network Analysis (WGCNA) [31]. The identified modules were characterized using GO-Elite [32] for controlling the network-wide false discovery rate. It was significant those pathways that had >10 genes at Z-summary >2 and FDR < 0.01 [31]. To estimate the enrichment of the DEG sets among the 17 modules, the Modular Single-set Enrichment Test (MSET) was applied with 10,000 permutations [33]. For that, we calculated the probability of the DEG being significantly represented in each module compared to 10,000 simulated sets of genes generated randomly from the microarray background. For each DEG and module, the p-value was calculated based on the number of simulated randomized sets that contains equal or higher number of genes belonging to the DEG in that module. Scatter smooth plots were construct only for the modules that had DEG over-representation (p-value ≤0.05) using the module eigengene values identified by the study [31].
In this study, the preservation module function from the WGCNA package in R were used to build unsigned co-expression networks to characterize features of coexpressed gene modules [34][35][36]. A pairwise correlation matrix was computed for all gene pairs across all samples, and an adjacency matrix was calculated by raising the correlation matrix to the power of 7. This threshold was chosen based on the scale-free topology criterion, which is the smallest threshold that resulted in a scale-free R 2 fit about 0.8, and used for all networks: schizophrenia samples only and controls samples only. For each pair of genes, a robust measure of network interconnectedness (topological overlap measure) was calculated based on the adjacency matrix. The topological overlap based dissimilarity was then used as input for average linkage hierarchical clustering. Finally, modules were defined as branches of the resulting clustering tree. To cut the branches, we used the hybrid dynamic treecutting in order to have robustly defined modules. The minimum module size was set to 20 genes with the deepSplit parameter set to 2. Each module was summarized by the first principal component of the scaled (standardized) module expression profiles (ME). Hub genes were assigned as having high intramodular connectivity kME >0.75. To discriminate changes in modules comparing SZP and CTS, we performed a module Preservation analysis that identified which modules were preserved in both datasets. For that, a Z-summary was calculated to exclude the possibility of randomness in preservation based on 1000 permutations, and indicating if a module is strongly (Z-summary >10), moderately (2 < Zsummary <10) or not preserved (Z-summary <2). For visualization, the resulting modules were plotted using Cytoscape [37].
For DEG, functional and disease enrichment analyses were performed in WebGestalt [38] using the whole genome as background with a further analysis implemented in DAVID (database http://david.abcc.ncifcrf.gov/). It was considered those biological processes and pathways that comprised more than 10 genes and FDR < 0.05.

Characterization of molecular alterations potentially interfering with gene expression
To refine the identification of genes regulated during neuronal differentiation, copy number alterations were characterized by array comparative genomic hybridization (CGH). Most of the detected copy number variations (CNVs) were present in all cell types (fibroblasts, hiPSCs and NPCs), which precludes their characterization as an acquired copy number alteration related to cell differentiation. There was no CNVs found exclusively in NPC from SCZP. Nevertheless, one large gain at 20q11.21 (1.2 MB) was acquired by the hiPSC and maintained in the NPC from the SCZP. The region contains 32 genes (Additional file 1: Figure S1). There were no CNVs found exclusively in NPC from SCZP.
To explore methylation alterations potentially involved in neuronal differentiation in schizophrenia, we compared NPCs, hiPSC and fibroblast from SCZP and CON, selecting CpG sites differentially methylated exclusively between hiPSCs and NPCs in SCZP and CON. A Venndiagram identified 51 (33 genes) and 5 (3 genes) CpG sites as hypo-or hyperrmethylated, exclusively in SCZP, respectively (Additional file 1: Figure S2, Table S1). These genes were over-represented in specific locations of the genome, 1p36, 10q22, 12q22 and 14q24 (adjusted p-value, adjP < 0.05). We also found overrepresentation of chromatin remodeling complex (adjP = 0.04), represented by DPF3, SMARCA2 and RERE. DPF3b interacts with the histones H3 and H4 in an acetylation-sensitive manner to regulate transcription [39]. SMARCA2 is a member of the SWI/SNF family of proteins with helicase and ATPase activities and thought to activate transcription altering chromatin structure in specific sites [40]. RERE is a nuclear receptor co-repressor that binds to the histone methyltransferase G9a orchestrating molecular events that lead to a stable methylation of histone H3-lysine 9 [41].

Characterization of expression patterns in SCZP
Gene expression profiles from hiPSC (SCZP or CON) were compared to the corresponding groups of differentiated cells (i.e. NPCs) to identify differences during neuronal differentiation. We considered that the analyses in CONderived cells identified DEG during neuronal differentiation; there were 1658 up-regulated (range from 0.50 to 6.76) and 1698 down-regulated (ranging from −0.43 to −8.49) genes. Comparison of the SCZP-derived cells revealed 752 up-regulated (range from 0.56 to 3.4) and 818 down-regulated (ranging from −0.54 to −4.97) genes. Subsequently, to exclude genes associated with normal neuronal differentiation, up-and down-regulated genes in CON and SCZP were compared to select those that were exclusive for SCZP. The analysis revealed 228 genes (151 up-and 77 down-regulated) altered only during neuronaldifferentiation of the SCZP cells (Additional file 1: Table S2, Figure S3), from which ten were described by studies in schizophrenia. Most of them were related to nervous system development and/or plasticity: NDEL1, ATF4, AKT1, MDGA1, PDE4B, FEZ1, TF3, ADRBK2, CACNA1C, ADRA2A (Additional file 1: Table S3). In agreement with our data, 10 genes (ATF4, CACNA1C, CPE, CTSF, EPHA2, GPR155, HOXB9, MYBL1, NLGN1, and PMAIP1) were identified by an independent study that used a similar design (hiPSC x NPC in schizophrenia cases) [26].
Using expression data from the prefrontal cortex of Brain Span (Allen Brain Atlas), expression of the 228 gene set were compared between 16 foetal and 22 nonfoetal samples. Nine genes were not present in Allen Atlas and 35 did not pass quality control (see material and methods). From the remaining 184 genes, 73 (40 %) were differentially expressed and therefore considered important for brain development.
To explore the functional significance of the genes from neurodevelopment associated with schizophrenia, we performed unbiased Gene Ontology term enrichment and pathway analyses of the 228 genes. The most noteworthy biological processes that were overrepresented were metabolic processes, signal transduction, nervous system development, regulation of neurogenesis and neuronal differentiation. Among the cellular pathways, we found over-representation of genes related to MAPK signaling pathway, pathways in cancer and metabolic pathways (Table 1). Finally, it was interesting to note that genes involved with cancer, infectious, metabolic and brain disorders (including stress, bipolar disorder and schizophrenia) were overrepresented diseases ( Table 2).
There were no overlap between the 228 DEG and genes from the 20q11.1 gain observed in the schizophrenic case. BACE1 was the only gene that was also identified by the methylation analysis. This gene seems to regulate proteins with important roles in developmental process such as NGR1 [42], which overexpression results in a schizophrenia-like behavior in a mouse model [43].
To gain further information of the 228 and 233 DEG sets, they were mapped to the gene co-expression modules from human cortical development [31]. The enrichment analysis point out that each set of DEGs had overrepresented genes but in different modules (M). While the set of 228 genes was enriched in M8, the set of 233 genes was enriched in M1, M5 and M15 (Additional file 1: Table S5). M8 and M15 were representative of genes preferentially expressed in "early" and "late" foetal period in human brain development, respectively (Additional file 1: Figure S4). This finding suggests that both datasets are related to neurodevelopment genes, but to different stages of the brain development.
To expand our findings from a cellular model of neuronal differentiation to postmortem samples from human brains, the 228 gene set, found exclusively in SCZP were compared to the 233 differentially expressed genes from adult brain samples. There were 194 genes out of the 228 represented in both platforms, from which four genes were common to both analyses: AQP1, DUSP6, FOS and TRIM24. Only DUSP6 and FOS have been described in schizophrenia [44,45]. None of the them were affected by the pH difference.
Further, to identify functional modules in patients and controls, a co-expression analyses of the 194 genes was performed in WGCNA. Modules for CTS and SZP were generated using the Scale-free Topology Criterion with power 7. Four modules were identified plus the grey module (Fig. 1a). Genes and biological processes overrepresented in each module are displayed in Tables 3  and 4, respectively. The best preserved module in SZP compared to CTS is the blue (blue, Z-summary >10, Fig. 1B) which had an overrepresentation of genes related to cell cycle.
Although showing a moderate Z-summary, using the median rank value based on the observed values (less dependent on modules sizes), the turquoise module was only marginally preserved (2 < Z-summary <10). Since this module includes negative regulation of cell differentiation and response to oxidative stress processes, identified deregulated in these same cells [24], Pearson correlation of paired genes were plotted (Fig. 1C). We can observe a change of correlations in SZP compared to CTS (change of color and thickness of lines linking two genes). Whilst co-expressed genes related to regulation of cell cycle processes, protein phosphatase binding and nonsense mRNA decay were extremely conserved between SZP and CTS, genes related to response to negative regulation of cell differentiation and oxidative stress were not.
Next, we defined the hub genes of the modules as those genes with high module connectivity (kME >0.75). This resulted in 13 hubs (CAPRIN1, CLIP4, CTNND1, GRB10, KHDRBS1, KPNA4, LYPD6, NRL, PSTK, RANBP2, RINT1, WASF3 and ZNF99) in the highly preserved blue module of CTS, all conserved in SZP (which gained additional four hubs: ZMIZ1, ORMDL1, EIF5 and FAM133B). It is interesting to note that in the yellow and brown modules, the number of hubs in SZP is smaller than in CTS, but all hubs identified in SZP were conserved from CTS, with only one exception in the yellow module (AKT) and one in the brown module (ZMPSTE24).
The marginally preserved turquoise module had 13 and 16 hub genes in CTS and SZP, respectively, with an overlap of only eight hubs, from which three were previously associated with schizophrenia: ADARB1 [46], PNOC [47] and XK [48]. Furthermore, three out of seven hubs found only in SZP were previously associated with schizophrenia. ADRA2A is required for normal presynaptic control of transmitter release in noradrenergic neurons and has polymorphisms associated to the increased risk of developing the disease [49]. CIT is a serine/threonine protein kinase involved in central nervous system development and cancer cells proliferation. Polymorphisms in CIT also seems to have an impact in the risk of developing schizophrenia [50]. LIN7A encodes a protein involved with the synaptic function and neuron migration in cerebral cortex development [51] and was associated to autism and attention deficit hyperactivity disorder [52,53]. Even though the other four hubs (ARMC8, BCL11A, CDO1 and FOXG1) have not been linked to schizophrenia, they are candidates for further studies.  Fig. 1 Module preservation analysis. a Dendrograms produced by average hierarchical clustering using topological overlapping matrix dissimilarity. Colours represent different modules. Upper panel shows modules in control subjects (CTS; blue, brown, yellow and turquoise) compared to patients with schizophrenia (SZP). Lower panel shows modules in patients with schizophrenia (blue, brown and turquoise) compared to modules in control subjects. b Left panel shows the median preservation rank (y-axis) in relation to module size (x-axis). Each circle represents a module labelled in different colours (blue, turquoise, yellow and brown). Right panel shows the Zsummary (y-axis) in function of module size. Dashed lines represent thresholds 2 and 10: ≥10: high preservation; 2 < Zsummary <10: moderate preservation; <2: low preservation. The panels show that the blue module is more preserved in control and patients whilst the turquoise module is less preserved. c Connectivity patterns (correlation network adjacencies) between genes from the turquoise module in controls (control) and patients (schizophrenia) showing a large loss of connectivity among genes in patient's module compared to control modules. Line thickness represents the connectivity pattern and line colour reflects the absolute correlation: −1 (negative) to 1 (positive). This module is enriched in genes related to response to negative regulation of cell differentiation and oxidative stress

Discussion
In the present study, we characterized molecular alterations during the differentiation of hiPSCs to NPCs from a patient with schizophrenia and health control and analyzed the identified genes in the adult brain tissue from schizophrenia cases and controls to study biological processes that could be disrupted in the disease. As we faced the limitation of a direct comparison due to the different nature of samples (single cell type-neuronversus all cells types-frontal cortex), we decided to make the analyses using gene co-expression network complemented with the exploration of biological processes. This allowed us to explore genes expressed during neuronal differentiation and important for the development of schizophrenia in the context of adult brain bulk.
It is known that cultured cells are more likely to have genomic instability, genetic and epigenetic mutations [54,55]. However, we only found one alteration by aCGH, a 1.2 MB gain in SCZP-hiPSC that was maintained in the NPC, suggesting that genes in this region may not be essential for the initial stages of neuron differentiation. Even though none of the 228 DEG was located in this region, it might contains regulatory variants with some role in the disease, although they were not explored in this study. The differentially methylated sites found during neuronal differentiation were preferentially associated with chromatin remodeling complex biological function, an interesting finding considering the previous study where the oxidative stress present in SCZP-NPC was reverted to NPC control levels after treatment of valproic acid (VPA), a histone deacetylase (HDAC) inhibitor [24].
Despite the fact that total RNA was extracted from three independent cell cultures (hiPSC and NPC), we are aware of the limitation of comparing one case with one control, which is the reason why we used the HT-self model for gene expression analyses. This consists in a mathematical model developed to analyze low numbers of samples [28]. Moreover, we used very stringent criteria in all analyses to identify the DEG and overrepresented ontologies (e.g. to be considered overrepresented, each class should contain more than 10 genes and a significant p-value after correction for multiple comparisons). Yet, data interpretations present certain limitations.
Concerning the genes involved with neuronal differentiation (hiPSC versus NPC), our comparisons revealed that cells derived from health control displayed twice more DEG than cells derived from the patient with schizophrenia. This finding may be involved with the disruption of cell differentiation that has been described by other studies in schizophrenia. The 228 gene set associated with neuronal differentiation in a patient with schizophrenia had over-representation of pathways classically related to cancer, such as MAPK and metabolic cellular pathways. Whilst new categories may appear with the advance of research in the field, many aspects of normal development and cancer have association to neuronal disorders [56]. There are no implications for tumorigenesis risk, although patients with schizophrenia are at a reduced risk of developing cancer [57]. In our previous study with these cells, we demonstrated that SZCP-NPC had altered metabolism compared to CON-NPC [24]. Likewise, the 228 set had over-representation of genes related to energy, lipid and carbohydrate metabolism, Table 3 Genes preserved in the modules of the respective samples Module Table 4 Biological processes overrepresented in the modules Module GO ID Ontology name Ontology definition Blue GO:0010564 regulation of cell cycle process Any process that modulates a cellular process that is involved in the progression of biochemical and morphological phases and events that occur in a cell during successive cell replication or nuclear replication events.  GO:0035264 multicellular organism growth The increase in size or mass of an entire multicellular organism, as opposed to cell growth.
GO:0042035 regulation of cytokine biosynthetic process Any process that modulates the frequency, rate or extent of the chemical reactions and pathways resulting in the formation of cytokines.
GO:0042089 cytokine biosynthetic process The chemical reactions and pathways resulting in the formation of cytokines, any of a group of proteins that function to control the survival, growth and differentiation of tissues and cells, and which have autocrine and paracrine activity.
GO:0042107 cytokine metabolic process The chemical reactions and pathways involving cytokines, any of a group of proteins or glycoproteins that function to control the survival, growth and differentiation of tissues and cells, and which have autocrine and paracrine activity.
GO:0006399 tRNA metabolic process The chemical reactions and pathways involving tRNA, transfer RNA, a class of relatively small RNA molecules responsible for mediating the insertion of amino acids into the sequence of nascent polypeptide chains during protein synthesis. Transfer RNA is characterized by the presence of many unusual minor bases, the function of which has not been completely established. which are altered in cancer [58] and stem cells [59] adding evidence for the implication of tumor-suppressor and oncogenes to neurodevelopment and associated disorders. Additionally, gene expression profiling from postmortem brain tissue of patients with schizophrenia and controls were studied. To begin with, although there were no significant differences in age or PMI, the pH was slightly lower in the brain of patients with schizophrenia, a characteristic that has been described although not considered for the analyses [15]. The pH may have an effect on the expression of some genes or these genes might be regulating the pH change. We found that 78 out of 233 DEG identified by the comparison between brain bulks from patients with schizophrenia and health controls had their expression also associated to differences in pH, an issue that needs to be deeper investigated.
Likewise, to explore the relationship between genes involved in neuronal differentiation and fully differentiated brains in patients with schizophrenia, we have to consider some additional known limitations: (1) Use of medication, hypoxia and/or injury are known confounding factors in the samples derived from the patients; (2) Comparison between data generated by single cell types (hiPSC and NPC) with data representing a mixed cell population at the end stage of the disease (a "whole brain tissue" that includes the glia, immunological and vascular cells); 3) Different methods were used to identify differentially expressed genes; (4) Genes from neurodevelopment are not being expressed or are being expressed at very low levels once the brain has ended its maturation. Supporting this later hypothesis, the comparison of the 228 and 233 DEG with modules derived from developing brains [31], showed that even though both sets have genes associated with neurodevelopment, they are representative of different stages. Thus, the overlap of only four genes was not a surprise.
Consequently, we decided to study the 228 gene set using co-expression modules in the expression data from the cohort of post-mortem cases of schizophrenia and controls identifying a well preserved module (blue) and, although not completely, a less conserved module (turquoise). The blue module was associated with regulation of cell cycle progress, protein phosphatase binding and nonsense mRNA decay. Genes from this module have function related to genome maintenance, RNA integrity, cell proliferation and/or differentiation, such as ANPC4, SPHK1 and SPP [60][61][62].
The turquoise module displays a set of genes involved in negative regulation of cell differentiation and oxidative stress response, more specifically reactive oxygen species (ROS) production, mitochondria function and neuron differentiation, which are closely related to the functional study of these cells that exhibited deregulation of mitochondria function in SZCP-NPC [24]. Among the genes represented in the module, three genes have a described function associated with ROS generation. ASK1 is a member of the mitogen-activated protein kinase family activated by cellular stress, activating c-jun N-terminal kinase (JNK) and p38 in response to oxidative stress, endoplasmic reticulum stress, infection and calcium influx [63]. PNPT1 is predominantly localized in the mitochondrial inter-membrane space and is implicated in RNA targeting to human mitochondria, small noncoding nuclear RNAs (5S rRNA, MRP RNA, some tRNAs, and miRNAs) [64]. cAMP-induced cFos is a transcriptional repressor of CD44 that triggers ROS generation [65]. Additionally, Table 4 Biological processes overrepresented in the modules (Continued) GO:0002764 immune response-regulating signaling pathway The cascade of processes by which a signal interacts with a receptor, causing a change in the level or activity of a second messenger or other downstream target, and ultimately leading to the activation, perpetuation, or inhibition of an immune response.
GO:0050865 regulation of cell activation Any process that modulates the frequency, rate or extent of cell activation, the change in the morphology or behavior of a cell resulting from exposure to an activating factor such as a cellular or soluble ligand.
GO:0002682 regulation of immune system process Any process that modulates the frequency, rate, or extent of an immune system process.
GO:0002684 positive regulation of immune system process Any process that activates or increases the frequency, rate, or extent of an immune system process.
GO:0010741 negative regulation of intracellular protein kinase cascade Any process that decreases the rate, frequency or extent of a series of reactions, mediated by protein kinases, which occurs as a result of a single trigger reaction or compound.
GO:0001932 regulation of protein phosphorylation Any process that modulates the frequency, rate or extent of addition of phosphate groups into an amino acid in a protein.
GO:0010627 regulation of intracellular protein kinase cascade Any process that modulates the rate, frequency or extent of a series of reactions, mediated by protein kinases, which occurs as a result of a single trigger reaction or compound.
some hub genes found exclusively in SZP have been associated to neuronal plasticity and different cancers, most of them associated to cell proliferation. Likewise, some of the hub genes have already been associated to neurodevelopment disorders. Mitochondria activity and energy metabolism undergo dramatic remodeling during embryonic development. Early embryos are initially dependent upon oxidative metabolism due to inheritance of maternal mitochondria which are subsequently segregated among daughter cells. A predominantly glycolytic metabolism provides sufficient energy to support stem cells in basal state, but a robust metabolic network is required to adequate the increasingly energetic demands [66]. The established initial mechanism support flux through the tricarboxylic acid cycle and electron transport chain whilst mitochondrial oxidization of pyruvate and glutamine to CO 2 is optimized to extract maximal energetic supply. It has been suggested that a concomitant rise in mitochondrial reactive oxygen species may prime stem cells for lineage differentiation [66,67]. The large number of observations identifying stem cell properties affected by energyresponsive molecules and signaling pathways raises questions about the fate of stem cells under conditions when metabolic homeostasis is perturbed. The literature suggests that abnormal endocrine signaling in organisms in extreme metabolic states has a substantial impact on proliferation and differentiation of multiple stem cell populations throughout the body [68]. It is also possible that variations in metabolism during gestation contribute to changes in offspring through their effects on stem cells. For example, low energy levels and the associated hormonal signals that occur in the pregnant mother could be directly transmitted to the offspring through the placenta, resulting in either transient or permanent changes in embryonic stem cells. Therefore, it seems that risk factors for schizophrenia such as famine, hipoxia and low birth weight could be implicated in permanent changes to stem cells leading to altered brain neural-development, as well as some alterations that are still presented in adulthoods, like disturbed neural plasticity and increased ROS production.