Variance of Gene Expression Identifies Altered Network Constraints in Neurological Disease

Using gene expression data from hONS and skin fibroblasts isolated from control donors, we next isolated the core attractor pathways whose differential expression distinguish a normal hONS stem cell expression phenotype, using the attract method that we recently developed [18]. Rather than testing individual genes, attract begins by using an ANOVA based method to test gene sets defined by KEGG pathways for their ability to distinguish between phenotypic states. Pathways ranked as significant are then each decomposed into “synexpression groups”—subsets with expression profiles that are both highly correlated and informative for distinguishing between phenotypes. These synexpression groups are then expanded to include genes with highly correlated profiles from within the original dataset, producing a collection of “core pathway modules.”

The top five attract modules were the MAPK signaling pathway, the focal adhesion pathway, the purine metabolism pathway, the neurotrophin signaling pathway and the cell cycle pathway (Table 1), each of which has previously been linked to important aspects of stem cell biology [19]–[21]. These modules were significantly different between the cell types tested at the 0.05 level (after adjusting for multiple testing using the Benjamini-Hochberg method). They were also the most “representative pathways,” in the sense that they contained the largest numbers of genes in the list of array probes significant between classes.

Overall we found 21 significant pathways (adjusted P-values <0.05). However, many of these pathways were overlapping in their gene content (see Table S1) and together represent three key common biological themes—immune response, growth factor signaling, and DNA replication—that are consistent with the phenotypic differences between the cell types.

In overlaying variance onto these networks we noted that the numbers of high-variance or highly constrained genes did not follow the expected population patterns, with a trend towards more high-variance genes across most of the pathways (Figure 3; Fig. 7 in Text S1). The purine metabolism and cell cycle pathways in particular had fewer low-variance and more high-variance genes than expected (Chi-square goodness of fit P-value <0.01 and P-value <1.2×10⁻⁶ respectively), an observation consistent in the hONS and skin fibroblast datasets. In contrast, the neurotrophin signaling pathway contained far fewer high-variance genes than expected (Chi-square goodness of fit P-value <0.01) suggesting that this network was under greater regulatory constraint.

Any analysis of variance must consider the contribution of technical variation to the data. To address this issue we compared the magnitude of the intra-individual expression variance with the size of the inter-individual expression variance associated with each donor group. Four additional donors had each contributed two independently derived biopsies, resulting in four replicate samples per donor (see Material and Methods and Text S3). The intra-individual mean of the CV distribution was smaller than the inter-individual CV mean, confirming that the technical variation in this data set was less than the biological variation observed. We next examined whether this observation held up when controlling for differences in sample size, by constructing variance distributions based on the same number of samples (n = 4), and found this trend persisted. Amongst the four individuals, differences in the CV distributions, albeit slight were observed for the five core pathway modules. Collectively these results suggest that repeated sampling of the same individual is associated with less expression variance, pointing to a strong genetic component in biological variability.

Disease status alters the variability profiles in core pathway modules

The distributions of expression variance for each of the patient groups highlighted an unexpected observation. Using a two sample t-test on log2 transformed CV values, we saw significant deviations from the control pattern in the SZ group (P-value <2.2×10⁻¹⁶), and also the PD group (P-value <0.001). This suggested that the SZ hONS lines demonstrated much less variation in their genome-wide expression patterns than was expected and in contrast, the PD hONS lines showed greater variation in their genome-wide expression patterns (Figure 4A; Text S2). Applying the same tests to identify disease-specific differences in genome-wide average expression showed no significant differences in either donor group.

We then investigated whether the differences in expression variance observed at a whole genome level were also apparent in our five core stem cell networks (Table 1). For each of these, we found that the SZ and PD groups sit at opposite ends of the variability spectrum (Figure 4C). The SZ group had a marked reduction in variance signaling pathway, as evidenced by more highly-constrained genes whereas the PD group had greater variance than the control group. The deviation in frequency distributions between the disease group and the control group was statistically significant for both SZ and PD groups in hONS stem cells, as assessed by a Chi-squared test with two degrees of freedom.

One might anticipate that increased expression variance, as seen in the PD group, is evidence of poor network integrity, a model that has been suggested as important in complex diseases. What was surprising, however, is the observation that the SZ group has significantly reduced expression variability in these core pathways, suggesting that both extremes of the variability spectrum may be implicated in disease processes.

These disease-specific patterns were independent on the particular cut-off used to define the regions of high and low constraint in the expression variance distribution (Figure 4). When more stringent cut-offs were applied (e.g. 5% and 10%) the same general trend of reduced variance for the SZ group and increased variance for the PD group was still observed.

Highly-constrained and lowly-constrained genes have different functional roles as reflected by their position in cellular networks

Disease status was associated with different proportions of lowly-constrained and highly-constrained genes; this observation raised the possibility that these two classes may play distinct roles in maintaining or driving cellular phenotype. Using Gene Ontology terms (cellular component, molecular function and biological process categories) [22], we performed a representation analysis for each set of highly-constrained and lowly-constrained genes in the pathway modules for the three donor groups. We then mapped these genes to networks based on highly-annotated protein-protein interaction data and compared the patterns of transcriptional constraint between phenotypic states.

We found that genes with high-variance/low constraints were both functionally and physically involved at the periphery of signal transduction pathways. In hONS cells, they functioned largely as cell surface receptors and tended to be localized in the membrane, transmembrane or extracellular matrix regions (Table S2). This might suggest that the hONS were heterogeneous for expression of growth factor receptors and as a population had dynamic interactions with the extracellular environment. In contrast, highly-constrained genes tend to function in signaling roles, such as protein kinases and phosphatases (see Text S4). This might imply that all of the cells in the hONS population were competent to transduce signals through the MAPK pathway, and were only restricted by the expression of receptors or the availability of ligands.

The hONS stem cells derived from SZ donors demonstrated both loss of variability at the cell surface, and increased constraint in the intracellular signaling molecules. For MAPK signaling, we saw significant functional enrichment of signaling GO categories (Table S2) for low-variance genes in the SZ and control groups, while there was enrichment for high-variance genes in the PD group. This suggests that MAPK signaling is particularly important for distinguishing between these three groups.

Just as cellular distribution of highly-constrained genes was not random, we observed a nonrandom pattern in the degree-distribution (connectivity) of genes based on their expression variance (Figure 5). The lowly-constrained genes had on average, a low degree whereas the highly-constrained genes were more highly connected to other genes in the network; this shift from a random distribution was marginally statistically significant for the control group (Chi-squared test; P-value  = 0.05958 for control group, P-value  = 0.1219 for SZ group, P-value  = 0.09595 for PD group, see Figs. 19 and 20 in Text S1), suggesting that not only are constraints imposed on genes linked to specific functional roles but they also have significantly distinct network topologies. Each disease group was associated with distinct deviations in degree distributions from those observed for the control group; the degree distributions between high and low-variance gene sets became more identical for the SZ group whereas in the PD group, the high and low-variance distributions appeared reversed from those observed in the control group. These observations suggest two different vehicles in which normal regulatory control through a signaling network may be disrupted or perturbed.

In comparing our three patient groups, we observed a consistent correlation between the degree of a particular gene and its expression variance, suggesting that transcriptional variability is an inherent aspect of all cellular systems (Fig. 20 in Text S1). This is consistent with our intuition that genes that are highly connected, and therefore play a central role in signaling or other networks, must be more tightly regulated than those that play more peripheral roles, including as cell surface receptors and downstream effectors.

It may be that the increased expressivity of those genes peripheral to the network provides a wider range of potential phenotypic response to external stimuli. In this way, disease processes that alter the variability in expression of particular genes may influence phenotype so that reduced expressivity limits the spectrum of response while increased expressivity tends toward loss of regulation over the cell's end-stage response.

Discussion

In genetics, the concepts of expressivity and penetrance describe phenotype variability between individuals with shared genotypes and across populations, respectively. The implicit explanations for phenotypic variability are differences in genetic, epistatic, or epigenetic interactions. In genomic analysis, we often see this variability in terms of sequence or structural polymorphisms within the genome [23] and focus our efforts on understanding the link between genetic variation and phenotypic diversity. In this context, low genetic variation leads to poor evolutionary fitness, whereas the convergence of multiple variants on disease networks is increasingly thought to contribute to disease states [24], [25].

A key element in the central dogma of molecular biology is the role played by RNA as an intermediary between gene and protein (and ultimately with phenotype). In this light it is surprising that variability of expression levels has received so little attention. This may reflect the fact that in the analysis of gene expression data, variance is often associated with technical artifacts rather than being seen as an intrinsic property that reflects the normal range of phenotypic heterogeneity. If technical noise intrinsic to the platform was the major driver in variance, then one would expect these random effects to affect experimental samples equally. Indeed, as our goal was to do comparisons between phenotypic groups, we made every effort to minimize experimental noise, beginning with standardized protocols for sample collection and laboratory handling, quality control during RNA extraction, and labeling and hybridization together so as to avoid potential batch effects.

Given that we see specific biases in the variance profiles that correlate in a meaningful way with phenotype argues that this is not the case. Our results suggest that this variability is, in fact, a much more important measure of phenotype than previously recognized and that changes in the spectrum of expressivity may be fundamentally linked to the development of distinct phenotypic states.

The expressivity of individual genes such as the pluripotency factors Oct4 and Nanog has been highlighted by others as an important contributor to phenotypic robustness of a population of embryonic stem cells. Variance of gene expression in these cells is both predictable and essential – providing a dynamic range of pluripotency factors which is directly linked to the differentiation potential of individual cells within that population [26], [27]. If expression variance of one or two key regulators is an important modifier of phenotype, then it holds that this can be measured at a pathway level as well. Even across cells derived from different donors, the distribution of expression variance for a pathway is predictable, and deviations from this correlate with disease phenotypes. What are the likely sources of this ‘regulated noise’? The reduced variance observed between hONS lines derived from different biopsies from the same donor indicates a role for genetics. One possible scenario is that genetic polymorphism might impose a combinatorial impact on the expressivity of genes within a network, and the consequent alteration of dynamic range of that network outcome. Indeed, the degree of expression variance of a reporter construct was found to be a heritable trait in S. cerevisiae [28]. Epigenetic factors are also likely to play a key role in expressivity at individual loci, and others have shown stochastic, population-wide variance in epigenetic modification of key developmental loci. With increasing evidence that expression variance is an important phenotypic attribute of cell populations, it holds that variance profiles may also reflect abnormal genetic or epigenetic events contributing to disease phenotypes.

In this context, it is perhaps not surprising that the most profound shifts in expression variance were found in hONS cells isolated from SZ donors. SZ is a life-long psychotic disorder, with age of onset in males in early adulthood and later in females. SZ is considered a disease of neurodevelopment, based on epidemiological, histological and genetic evidence [29]–[32]. It is clear that SZ is a complex genetic disease with a strong environmental component. There are several well replicated genome-wide studies that have implicated common polymorphisms in the etiology of the disease, although these account for a small component of the heritability and an emerging theory is that polygenic risks explain more of the genetic component of this disease. Reduced expression variance of hONS networks provides fresh clues into potential mechanisms underlying diseases like SZ. For example, if the patterns of gene expressivity found in the hONS reflect those found in the developing brain, then the orderly cascade of brain development may be altered. In addition, it is feasible that an overly constrained biological pathway would be less adept at buffering environmental stressors.

Neuronal stem cells are the progenitors for neural cell states and the development of brain is rooted in the cell fate decisions that occur in these stem cells. Olfactory stem cells have been shown to serve as a good surrogate for neuronal stem cells and stem cell differentiation[14]. Our analysis of gene expression in olfactory stem cells cultured from patients with SZ, PD, and matched controls identified a number of key pathways that distinguished these three groups, including key signaling and developmental pathways. What was surprising about these pathways was that not only was there a notable difference in expression, but that expressivity, or variability in gene expression levels, was also significantly different between these different diseases. Further, when the expressivity was mapped to protein-protein interaction networks, there were distinctly different patterns of transcriptional constraint that depended on the connectivity of the proteins in the network. Further, these patterns depended not only on disease state, but also on the degree of connectivity within the corresponding protein-protein interaction pathway.

In examining the network topology, we found that in the SZ patients, there was significantly greater regulatory control over genes at the highly-connected core of the corresponding pathways. For example, in the MAPK pathway, SZ patients exhibited far less expressivity in the kinases that represent the core of the pathway. On the other hand, the PD patients had significantly fewer constrained genes mapping to the same core pathway regions than did controls. The increased variance associated with some proteins such as receptors, may indicate that these proteins have fluctuating turnover rates in the cell populations, which in turn will influence the capacity of cells to interact with their environment.

These data, together, paint a very interesting and compelling picture of the mechanism associated with disease. One of the defining characteristics of SZ is the disruption of normal cognitive processes and this is reflected in fewer neural synapses in parts of the brain. While it is no doubt a leap to interpret our results as having a direct link to thought processes in the disease, one could imagine that a high degree of constraint in key transcription factor networks that play a role in cellular differentiation and developmental processes. One could imagine that highly constrained transcriptional networks in neuronal stem cells could reduce the spectrum of cellular phenotypes that could be derived from those stem cells and reduce the plasticity resulting in the brain and altering its potential responses and thought patterns.

It is worth noting that among the pathways that demonstrate highly-constrained gene expression in our SZ patients are those involved in signaling in cancer. One consequence of increased constraints, and reduced plasticity, is defects in the network are more likely to be disastrous. In this kind of model, defects in cancer-related networks would lead to loss of those cells, preventing the types of adaption observed in cancer development. Although there are conflicting reports regarding the risk of cancer in schizophrenics, the majority of reports suggest that patients with SZ are protected against cancer in general, and from lung and colorectal, despite increased smoking [33], [34] and drinking habits [35] in this population. Laboratory studies have also reported reduced tumor growth in animal models of schizophrenia [36]. Although speculative, the high degree of transcriptional constraint we found in the MAPK pathway in SZ patients may explain, in part, these observations.

In contrast, disease states may arise where increased variance changes the predictability of network outcomes, resulting in dysregulation of the desired cell state. In PD, we observed an increase in the expression variance of core signaling pathways, which we predict will diminish the robustness of the network to external events. This may be an essential element that is shared amongst diseases of aging.

Although SZ and PD represent very different conditions, our results suggest that changes in expressivity relative to the normal spectrum of variability, may play an important role in the development of disease phenotypes. One possible way to interpret this is in the context of models first proposed by Conrad Waddington [37] and later refined by Stuart Kauffman [38]. Waddington and Kauffman envisioned what we can interpret as a gene expression state-space landscape defined by possible gene expression states. In this landscape, stable cell states represent fixed points (Kauffman referred to these as “attractors”) connected by evolutionarily-defined “canals” representing the differentiation pathways connecting distinct cell phenotypes. In this model, the SZ patients would be characterized by more highly-constrained, “deeper,” canals, limiting the potential end states that one might achieve during differentiation. PD patients, on the other hand, could be characterized by a flattening of the same canals (what has been referred to as “decanalization”), increasing dysregulation of the associated pathways and potentially allowing for a degradation of the well-defined cellular end states that might be available. While the potential importance of decanalization in disease has been discussed [25], we believe that our results are the first to suggest that over-canalization may be an equally important process in developing disease phenotypes.

Materials and Methods

Ethics statement

This work uses public expression datasets from patient-derived cells which were collected under the ethics approval of the Griffith University ethics committee.

Gene expression data set

Illumina Human-Refseq8 v2 BeadChips (Illumina, Inc.) arrays were used to capture genome-wide gene expression profiles; the raw data was summarized using BeadStudio Version 3.1.7(Illumina, Inc). Background correction and normalization methods were performed using the R/Bioconductor lumi package. All downstream analyses were performed using Quantile normalized data, without background correction, and only probes passing the Illumina detection threshold were included in variance analysis. A probes was considered to pass the Illumina detection score if it had a detection p-value ≤0.01 in at least 75% of individuals in the same donor group, resulting in 14,986 probes. This expression data is available from ArrayExpress under the experiment accession number E-TABM-724.

Intra-individual variance

Illumina Human-Refseq8 v3 BeadChips (Illumina, Inc.) arrays were used to collect genome-wide gene expression profiles of the four replicate samples for each of the four donors. The donors were made up of two healthy controls and two PD patients, and each donor underwent two independently derived biopsies that were each replicated twice, giving rise to four samples per donor. The raw gene expression data set was summarized as above. The data set was filtered using the same detection filter, and the subset of 6,809 detected probes common to the v3 and v2 arrays was used for the comparison of intra-individual and inter-individual variance analysis. To gauge the effect of sample size on variance, the inter-individual variance was calculated by computing CV distributions based on a random subset of 4 individuals from each donor group, and comparing these to the intra-individual CV distributions that were calculated from the four samples for each individual. From the 100 random subsets generated, we observed a reduction in the difference between the inter-individual variance and the intra-individual variance when the sample size was reduced and fixed at four replicates.

Attract method

The attract package can be obtained from Bioconductor [http://www.bioconductor.org/packages/devel/bioc/html/attract.html] and is available as a module in the MeV microarray analysis software [39] (http://www.tm4.org/mev). To identify activated core pathways whose expression defines a control hONS phenotype, we ran attract on an expression data set consisting of the skin fibroblasts and two types of hONS lines from the control patients and a set of mesenchymal stem cell lines from a group of unrelated individuals. Attract was run using pathway modules defined by KEGG pathways represented in Bioconductor (version 2.4.1, Biobase version 2.8.0 and illuminaHumanv2BeadID.db version 1.6.0).

Identifying low- and high-expression variance genes

A CV value was calculated for each detected probe by dividing the standard deviation of its expression in a donor group by its average group expression. Low and high-expression variance genes were identified as those genes below and above the 25^th percentile of the genome-wide CV distribution based on values from all donors.

Assessing the significance of disease status versus healthy control group

P-values were obtained by comparing the counts of high, medium and low constraint genes in each of the control, SZ and PD groups and using a Chi-squared goodness of fit test where counts from the control group were designated as the expected counts.

Representational analysis

For each of the top five pathway modules, we applied a representational analysis to the set of highly-constrained and lowly-constrained genes for each of the donor groups. We used tools from the Bioconductor GOstats package (version 2.14.0, run on R version 2.11.1). P-values were adjusted for using the Benjamini-Hochberg method within each ontology class, donor group and pathway module (Cellular Component, Molecular Function, or Biological Process) and significant results obtained at the 0.05 level. To focus on what functional terms were unique to genes of altered constraint, we excluded significant GO terms that appeared in both lists of highly-constrained and lowly-constrained genes and retained only those GO terms that were unique to each list. The significant GO terms for highly-constrained genes appear in the Supplementary File.

Network topology p-values

A Kolmogorov-Smirnov test was used to test the degree distributions of the lowly and highly-constrained gene sets. A Gaussian kernel density estimator (the density function from R, using the default method to select bandwidth size) was used to produce the degree distribution density plots shown in Figure 4.

Literature-curated networks

Protein-protein interactions were defined using two knowledge-based annotation systems, Ingenuity Pathway Analysis (IPA) software and the GeneGo metacore tool. Both tools permitted identification of highly curated protein-protein interactions; using IPA we extracted the degree of connectivity for each gene in the attract-networks and an image of the interaction network was obtained from GeneGo.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6