Parkinson's disease and schizophrenia interactomes contain temporally distinct gene clusters underlying comorbid mechanisms and unique disease processes

Abstract

Genome-wide association studies suggest significant overlaps in Parkinson's disease (PD) and schizophrenia (SZ) risks, but the underlying mechanisms remain elusive. The protein-protein interaction network ('interactome') plays a crucial role in PD and SZ and can incorporate their spatiotemporal specificities. Therefore, to study the linked biology of PD and SZ, we compiled PD- and SZ-associated genes from the DisGeNET database, and constructed their interactomes using BioGRID and HPRD. We examined the interactomes using clustering and enrichment analyses, in conjunction with the transcriptomic data of 26 brain regions spanning foetal stages to adulthood available in the BrainSpan Atlas. PD and SZ interactomes formed four gene clusters with distinct temporal identities (Disease Gene Networks or 'DGNs'1-4). DGN1 had unique SZ interactome genes highly expressed across developmental stages, corresponding to a neurodevelopmental SZ subtype. DGN2, containing unique SZ interactome genes expressed from early infancy to adulthood, correlated with an inflammation-driven SZ subtype and adult SZ risk. DGN3 contained unique PD interactome genes expressed in late infancy, early and late childhood, and adulthood, and involved in mitochondrial pathways. DGN4, containing prenatally-expressed genes common to both the interactomes, involved in stem cell pluripotency and overlapping with the interactome of 22q11 deletion syndrome (comorbid psychosis and Parkinsonism), potentially regulates neurodevelopmental mechanisms in PD-SZ comorbidity. Our findings suggest that disrupted neurodevelopment (regulated by DGN4) could expose risk windows in PD and SZ, later elevating disease risk through inflammation (DGN2). Alternatively, variant clustering in DGNs may produce disease subtypes, e.g., PD-SZ comorbidity with DGN4, and early/late-onset SZ with DGN1/DGN2.

Fig. 1. Flow chart depicting the methodology and the findings of the study.

To examine the neurobiological underpinnings of PD and SZ, the study focused on disease interactome analysis. Genes associated with PD and SZ were compiled and used to construct the PD and SZ disease interactomes. Next, the temporal expression profiles of the genes in these networks were examined using hierarchical clustering to identify temporally distinct clusters. Four temporally distinct clusters were extracted from the disease interactomes (disease gene networks or ‘DGNs’ 1-4). Their temporal characteristics are listed under “inferences from temporal expression profiles”. Their enrichments for specific sub-networks in the disease interactomes are listed under “enriched disease interactome sub-networks”. Finally, the four DGNs were characterized based on their enrichment patterns in specific brain regions, cell types, pathways, GWAS traits, and disease endophenotypes curated from literature. Their relationships with the interactomes of 22q11.2 deletion syndrome and Lewy body dementia, in which affected individuals develop both psychosis and Parkinsonian symptoms, were also examined. The overall findings from these analyses are listed under “inferences from functional enrichment analyses”.

Fig. 2. Genes in PD and SZ interactomes split into four clusters based on their temporal expression patterns.

a The figure shows the two spatiotemporal clusters – corresponding to 16 brain regions in early infancy to adulthood stages (left) and 25 brain regions in early prenatal to late prenatal stages (right) – on the horizontal axis. The relative expression of 4436 out of the 4629 genes present in the PD and SZ interactomes across 407 spatiotemporal points was hierarchically clustered. Clustering was performed on log-transformed RPKM (i.e. Reads Per Kilobase per Million mapped reads) values using the hierarchical clustering method with average linkage. The dendrograms were derived from the clustering analysis based on the computation of Pearson correlation coefficients between the data points. The clustered heat map was created using the Morpheus software. Two primary gene clusters can be seen on the vertical axis, which can be further subdivided into five sub-clusters (A-E) showing distinct temporal profiles. These sub-clusters were labelled as Disease Gene Networks or ‘DGNs’1-4 based on their preferential enrichment for proteins uniquely found in the PD/SZ interactomes or shared between both the interactomes. b–d PCA was performed with the expression profiles of the genes belonging to PD and SZ interactomes in 2 brain regions in various developmental stages (A1C-24pcw, A1C-25pcw, A1C-37pcw and AMY-21pcw). A matrix with 4436 genes (rows) and 4 spatiotemporal conditions (columns) was constructed out of log-transformed RPKM values. Unit variance scaling was applied across this matrix. Singular value decomposition (SVD) with imputation was used to extract principal components (PCs). Component scores corresponding to PC1’ and PC2’ explaining 89.5% and 6% of the total variance were plotted along X and Y axes respectively. Component scores and quadrant-wise enrichment of genes (b) shared between PD and SZ interactomes, (c) uniquely found in the PD interactome, and (d) uniquely found in the SZ interactome. The percentage of genes from each of these gene sets found in each of the 4 quadrants is also shown. Note that these analyses were performed using data collected from a single source (Human Developmental Transcriptome, BrainSpan Atlas).

Fig. 3. DGNs 1-4 showed differential enrichment patterns in specific brain regions, signalling pathways and brain cell types.

The figures show the enrichment of the DGNs for genes (a) expressed across 26 brain regions compiled from BrainSpan Atlas having logRPKM > 2, (b) showing high/medium expression in 13 brain regions (Transcripts Per Million or TPM > 9) compiled from GTEx, (c) involved in five neuronal synaptic signalling pathways, i.e., cholinergic, dopaminergic, GABAergic, glutamatergic and serotonergic pathways compiled from the KEGG database, and (d) specifically expressed in 6 brain cell types compiled from a study by Lake et al. The statistical significance of each enrichment was computed as a p-value. All p-values were transformed to –log₁₀p-values, and then assembled into a data matrix containing brain regions as rows and the DGNs as columns. Negative log transformation simplifies the scale, revealing patterns and significant results more visibly, with higher –log₁₀ values signifying smaller p-values. Variations in enrichment are represented in the form of heat maps. Specifically, z-scores computed based on the inverse normal transformation of –log₁₀ transformed p-values are shown. Clustering was performed using the hierarchical clustering method with average linkage. The dendrograms were derived from the clustering analysis based on the computation of Pearson correlation coefficients between the data points. The clustered heat map was created using the ClustVis software. Note that this analysis was performed with gene expression profiles of brain samples obtained from a single source (BrainSpan Atlas).

Fig. 4. Lewy body disease and 22q11 deletion syndrome interactomes clustered with specific DGNs.

The enrichment of the genes belonging to DGNs, and the interactomes of GWAS genes associated with Lewy body disease (LBD) and genes differentially expressed in SZ patients with 22q11 deletion (22q11 del), in (a) 26 brain regions from BrainSpan Atlas and (b) 13 brain regions in GTEx. The statistical significance of each region-wise enrichment was computed as a p-value. All p-values were transformed to –log₁₀p-values, and then assembled into a data matrix containing brain regions as rows and the DGNs as columns. Negative log transformation simplifies the scale, revealing patterns and significant results more visibly, with higher –log₁₀ values signifying smaller p-values. Variations in region-wise enrichment have been represented in the form of heat maps. Specifically, z-scores computed based on the inverse normal transformation of –log₁₀ transformed p-values are shown in (a, b). Clustering was performed using the hierarchical clustering method with average linkage. The dendrograms were derived from the clustering analysis based on the computation of Pearson correlation coefficients between the data points. The clustered heat map was created using the ClustVis software. Note that the results shown in (a, b) were obtained using data from a single source, i.e., BrainSpan Atlas in the case of (a) and GTEx in the case of (b).