Dynamic perturbations of cellular processes associated with prion diseases
Previously, we proposed core genes associated with pathogenesis of prion disease that showed shared dynamic changes in differential expression along progression of prion disease [2]. To identify these core genes, we first selected the following three combinations among the eight prion strain-mouse strain combinations previously reported [2]: 1) inbred C57BL/6 J (B6) and 2) FVB/NCr (FVB) mouse strains infected with RML (B6-RML and FVB-RML), and 3) FVB background mice expressing half of the amount of PrPC (Prnp0/wt) as wild type mice infected with RML (Prnp0/wt-RML). Prnp0/wt-RML mice have a longer incubation times (more than 250 days) than B6- and FVB-RML (~ 150 days) and a slower rate of synaptic degeneration than B6- and FVB-RML, despite a rate of PrPSc accumulation similar to that of FVB-RML. The inclusion of Prnp0/wt-RML enabled us to reliably discern early and late responsive genes.
Using the time-course gene expression profiles of the three selected prion-mouse strain combinations (B6-RML, FVB-RML, and Prnp0/wt-RML), we performed the orthogonal non-negative matrix factorization (ONMF) clustering [34] with the number of clusters = 20. Among the 20 clusters (Additional file 1: Figure S1), we focused on three clusters (Clusters 1–3) including 883 differentially expressed genes (DEGs) that showed the following shared differential expression patterns across the three combinations: 107 early up-regulated genes in Cluster 1 (Fig. 1a, EU); 502 late up-regulated genes in Cluster 2 (Fig. 1b, LU); and 274 down-regulated genes at the late stage in Cluster 3 (Fig. 1c, DN). Many of these 833 DEGs have been identified previously as DEGs in prion-infected mice (Additional file 2: Table S1) [35,36,37,38,39,40,41,42]. Moreover, of the 883 DEGs, 270 and 541 overlapped with the previous 333 core genes and 923 prion disease-related genes identified from five prion-mouse strain combinations (B6-RML, FVB-RML, C57BL/6.I-1(B6.I)-RML, B6-301V, and B6.I-301 V), respectively (Additional file 2: Table S1).
To examine how these genes represent early and late perturbations of cellular processes associated with prion diseases, we performed enrichment analyses of gene ontology biological processes (GOBPs) for EU, LU, and DN genes (Fig. 1d; Additional file 2: Table S2) using DAVID software [29]. The GOBPs represented by each cluster of genes reflect cellular processes perturbed with the corresponding dynamics. Early activated GOBPs (Groups 1 and 2 in Fig. 1d) include the processes related to innate immune or inflammatory responses (e.g. antigen processing and presentation, complement activation, pattern recognition receptor signaling pathway, phagocytosis, and leukocyte activation). Late activated GOBPs (Group 3 in Fig. 1d) include: 1) innate and/or adaptive immune responses (e.g. cytokine production, cytokine-mediated signaling, NFKB/MAPK signaling, and B-cell activation); 2) PrPSc deposition and transport (vacuole organization, regulation of proteolysis, extracellular matrix organization, and actin cytoskeleton organization); 3) lipid metabolism and transport (sphingolipid metabolism and cholesterol transport); and 4) stress responses (angiogenesis, responses to reactive oxygen species, and apoptosis). Finally, down-regulated GOBPs (Group 4 in Fig. 1d) include synaptic degeneration-related processes (axon guidance, calcium ion transport, synaptic transmission, and gamma-aminobutyric acid signaling). Early and late perturbations of these processes were consistent with those in the previous PPI network models [2].
TRN describing transcriptional regulation of early and late perturbed processes associated with prion diseases
To investigate transcriptional regulation underlying early and late perturbations of cellular processes associated with prion diseases, we built TRN describing regulations for EU, LU, and DN genes by 467 TFs (Additional file 1: Figure S2A) based on TF-target interactions in MetaCore™ (ver 6.7) and Ingenuity Pathway Analysis (IPA). Based on these TF-target interactions, the upstream TFs were available for 84.7% (748 of 883) of the EU, LU, DN genes. The 467 TFs in the TRN were first categorized into two groups based on their differential expression: 1) 40 differentially expressed TFs (DETFs) including 5 EU, 24 LU, and 11 DN TFs and 2) 427 non-differentially expressed TFs (non-DETFs). Of the 107 EU genes, 90 (84.1%) can be transcriptionally regulated by 26 DETFs and 160 non-DETFs (Fig. 2a). Of the 502 LU genes, 448 (89.2%) can be regulated by 36 DEFTs and by 368 non-DETFs (Fig. 2b). Finally, of the 274 DN genes, 210 (76.6%) can be regulated by 24 DETFs and 177 non-DETFs (Fig. 2c). Interestingly, the DETFs can regulate larger numbers of the target genes per TF compared to the non-DETFs (Figs. 2a-c). For example, 66 EU genes were regulated by 26 DETFs (66/26 = 2.54), while 85 EU genes were regulated by 160 non-DETFs (85/160 = 0.53) (Fig. 2a). These data indicate greater relative importance of the DETFs in transcriptional regulation of the EU, LU, and DN genes.
Among the 467 TFs, we next selected 114 major TFs that have significant (FDR < 0.1) numbers of target genes in EU, LU, or DN clusters (Additional file 2: Table S3). Based on the observation that DETFs have larger numbers of targets than non-DETFs, among the 114 major TFs, we selected the following 18 DETFs as key TFs that play important roles in transcriptional regulation of cellular processes represented by the EU, LU, and DN genes: 1) 4 EU-DETFs (Cebpa/d, Irf8, and Nupr1); 2) 12 LU-DETFs (Atf3, Cebpb, Jun, Nfe2l2, Rela, Spi1, Srebf1, Stat3/6, Ikzf1, Nfkbia, and Id3); and 3) 2 DN-DETFs (Atf2 and Sox2). Interestingly, these 18 key TFs, including two down-regulated TFs, were found to regulate mostly the EU and LU genes, according to the TF-target interactome data used (Fig. 2d). Several of the 18 key TFs have been reported to have potential involvement in prion diseases or other neurodegenerative diseases (Additional file 2: Table S3). For example, Stat3 showed increased phosphorylation in mouse brains infected with prions, suggesting its potential roles in pathogenesis of prion disease [43]. Also, Atf2, an abundant TF in normal brains, was significantly down-regulated in the brains with Alzheimer’s, Parkinson’s, and Huntington’s diseases, consistent with our findings [44]. Moreover, a dominant negative TF Jun significantly reduced neuronal death in prion infected neurons [45], and Nfe2l2, also known as Nrf2, was strongly up-regulated in multiple sclerosis lesions and found to be associated with active demyelination in the lesions [46].
We then developed a TRN model using the 18 key TFs and 315 target genes (42.1% of 748 DEGs) based on the TF-target interactome data (Fig. 2e). The target genes in the TRN were organized into three subnetworks that represent three aspects of disease: PrPSc accumulation, microglial/astrocytic activation, and synaptic degeneration [2]. We then examined how significantly the 18 key TFs regulate cellular processes associated with the three pathological features (Fig. 1d) by computing the fraction of their target genes to the DEGs involved in each cellular process (Additional file 1: Figure S2B; Additional file 2: Table S4). For example, the processes represented only by the EU genes (Group 1 in Additional file 1: Figure S2B) are regulated by 10 of the 18 TFs (Cebpa/d and Irf8 EU TFs; and Cebpb, Jun, Rela, Spi1, Stat3, Nfkbia, and Atf2 LU TFs). The processes represented by EU and/or LU genes (Groups 2 and 3 in Additional file 1: Figure S2B) are regulated by most of the 18 TFs. The processes represented by the DN genes (Group 4 in Additional file 1: Figure S2B) are regulated mainly by eight of the 18 TFs (Jun, Rela, Atf2, Cebpb, Stat3, Spi1, and Srebf1 LU TFs; and Sox2 DN TF). Taken together, the TRN delineated by the 18 key TFs and the 315 target genes effectively covered transcriptional regulation of early and late cellular processes associated with the pathological features of prion diseases (Fig. 2f).
Over-represented regulatory motifs in the prion disease-associated TRN
Several methods have been developed for analyzing regulatory structures of biological networks to identify core regulatory motifs or modules in the networks, including significant area search, network propagation, clustering-based methods, and network motif analysis [6]. To determine a method suitable for our TRN, we first examined the topological characteristics of the original TRN constructed for all 467 TFs by analyzing distributions of degrees and clustering coefficients for nodes. The distributions showed that the TRN has the features of both scale-free and hierarchical networks (Figs. S3A-B), indicating that the TRN may include functional regulatory motifs [33]. Collective operations of network motifs composed of TFs and their targets characterize early and late induction or repression of target genes. Thus, among the previous methods, we employed the analysis of network motifs (Fig. 3a, top left) to examine detailed transcriptional regulatory structures in our TRN associated with early and late alterations of target genes.
For network motif analysis, we performed the enrichment analysis of 13 different previously reported regulatory motifs, each of which included a pair of TFs and a target gene (TFi-TFj in Fig. 3a, upper left), as previously described [47]. In this analysis, we used two TRNs constructed for all 467 TFs in addition to the 114 major TFs, not for the TRN for the 18 key DETFs to avoid bias toward the key DETFs in the enrichment analysis. The analysis revealed that among the 13 regulatory motifs, only two motifs (Motifs 7 and 10) including feed-forward loops were significantly (P < 0.01) over-represented consistently in the two TRNs (Fig. 3a). Thus, we used the TRN for the 114 major TFs for the following analysis to focus on Motifs 7 and 10 defined by the major TFs. The TRN included 6329 Motif 7 and 533 Motif 10 that comprised 95 and 15 major TFs, respectively, and large portions of the DEGs as target genes (56 EU, 228 LU, and 54 DN genes for Motif 7; and 30 EU, 125 LU, and 23 DN genes for Motif 10) (Fig. 4b; Additional file 2: Table S5). Interestingly, Motifs 7 and 10 were also found to be over-represented in the TRNs for E. coli, yeast, fly, and/or sea urchin [47], suggesting that the TRN in prion disease includes key regulatory structures conserved in mammalian TRNs in spite of the incomplete nature of TF-target data.
We showed above that the DETFs have higher contributions to regulation of target genes than non-DETFs (Figs. 2a-c). Thus, among the TF pairs in Motifs 7 and 10, we next selected 48 and 20 DETF pairs from 6329 Motif 7 and 533 Motif 10, respectively (Fig. 3c and Additional file 1: Figure S3C, 2nd box). Finally, among them, we identified 20 and 15 key DETF pairs that were significantly (FDR < 0.1) enriched in the 6329 Motif 7 and 533 Motif 10, respectively (Fig. 3c and Additional file 1: Figure S3C, 4th box). Interestingly, these selected key DETF pairs (Fig. 3d) included the 18 key DETFs (Fig. 2d), targeting different sets of 178 genes in Motifs 7 and 10, respectively. To understand the contribution of the key DETF pairs to the regulation of cellular processes related to prion disease, we then examined how significantly their target genes in Motifs 7 and 10 overlapped with the DEGs involved in cellular processes (Fig. 1d) associated with the pathological features in prion diseases. A majority of the key DETF pairs showed significant (P < 0.01) overlaps of their targets with the DEGs involved in cellular processes related to PrPSc replication and accumulation and microglia and astrocyte activation, suggesting their strong associations with these pathological features (Fig. 3d and Additional file 1: Figure S3D).
Regulatory motifs that can be operative in individual cells of the brain
In diverse types of cells, such as astrocytes, microglia, oligodendrocytes, and neurons in the brain, molecular networks may undergo alterations during the progression of prion disease [7]. Our gene expression data were generated from the whole brain, providing mRNA expression levels from the mixture of diverse cells in the brain. Thus, it is difficult to sort out in what cell types Motifs 7 and 10 were altered by prion infection. To examine whether Motifs 7 and 10 can be operative in individual cell types, we first analyzed whether both TFs in each selected TF pair were expressed in the same cell type using previously reported gene expression data of seven cell types in the mouse brain, including microglia, astrocytes, neurons, oligodendrocyte precursor cells (OPCs), newly formed oligodendrocytes (NFOs), myelinating oligodendrocytes (MYOs), and endothelial cells [48]. For example, for Cebpa-Jun pair forming Motif 7, both Cebpa and Jun were found to be expressed only in microglia, neurons, and OPCs among the seven cell types (Fig. 4a). To evaluate whether TFs were expressed, FPKM ≥5 was used as a cutoff as in Zhang et al [48]
We then examined whether such expressed TF pairs had significant regulatory power by analyzing the numbers of their target genes expressed in the same cell types. The Cebpa-Jun pair had 51 targets in the 6329 Motif 7, and significant (FDR < 0.1) numbers of the target genes were found in all three cell types (14, 28, and 37 genes in neuron, OPCs, and microglia, respectively) (Fig. 4b), suggesting that Motif 7 formed by Cebpa-Jun pair and their targets can be operative with significant regulatory power in these three cell types. Finally, we analyzed the correlation between differential expression patterns of TFs and their targets, which were expressed in the same cell type, in the whole brain data for the B6-RML combination along disease progression (Fig. 4c). For the Cebpa-Jun pair, Cebpa and Jun first showed a significant (P < 0.01) positive correlation (0.84) between their differential expression patterns along the progression of prion disease (Fig. 4d). Additionally, the targets expressed in OPCs and microglia, yet not in neurons, showed significant positive correlations (≥ 0.6 statistically) with Cebpa and Jun (Fig. 4d). All these data suggest that Motif 7 formed by Cebpa-Jun and its targets can act as a functional motif with significant regulatory power in OPCs and microglia during the progression of prion disease.
We applied the same procedure to all 20 and 15 key DETF pairs selected for Motifs 7 and 10, respectively, and identified 8 and 6 key DETF pairs for Motifs 7 and 10, respectively, which could be operative with significant regulatory power in one of the seven cell types (Fig. 4e and Additional file 1: Figure S4). Interestingly, these DETF pairs were found to have significant regulatory power most strongly in microglia (9 and 3 pairs for Motifs 7 and 10, respectively) and OPCs (7 and 2 pairs for Motifs 7 and 10, respectively). Only few DETF pairs had significant regulatory power in astrocytes (3 pairs), neurons (1 pair), MYO (0 pair), and endothelial cells (3 pairs). Moreover, the comparison of the DETF pairs between microglia and OPCs revealed that the target gene counts for the DETF pairs were higher in microglia than in OPCs (Fig. 4b). Collectively, all these data suggest that these key regulatory motifs are employed in microglia and OPCs to regulate the target genes in the two cell types. Moreover, they are more strongly operative in microglia than in OPCs.
Core transcriptional regulatory circuits for dynamic activation of cellular processes
To understand collective actions of the regulatory motifs, we next combined all the regulatory Motifs 7 and 10 that could be operative in microglia and OPCs with strong regulatory power and identified core transcriptional regulatory circuits (TRCs) including key DETF pairs and their target genes in the selected Motifs 7 and 10 in the two cell types (Figs. 5a-b). In the microglial and OPC TRCs, EU and LU genes were targeted by EU or LU DETFs and found to be more regulated by both EU and LU DETFs than by either EU or LU DETFs alone (Figs. 5c-d). Interestingly, for the LU target genes, their fold-changes tended to be higher when regulated by both EU and LU DETFs than by either EU or LU DETFs alone (Figs. 5e-f). All these data suggest that EU DETFs modulate the expression of LU target genes, and LU DETFs also modulate the expression of EU target genes at the late stage of prion diseases.
Comparisons of the nodes in the microglial and OPC TRCs showed that 8 DETFs and 18 target genes were shared between the two TRCs, which corresponded to 80% of the DETFs and 43% of the target genes in the microglial TRC, suggesting that the same DETFs target different sets of genes between microglia and OPCs (Figs. 5a-b). The target genes in the microglial TRC were mainly involved in early inflammatory and immune responses (Fig. 1d) represented by the EU genes (Fig. 5g; Additional file 2: Table S6). Compared to the microglial TRC, the OPC TRC included the target genes that were mainly involved in the relatively later immune responses and responses to oxidative stress and steroid hormones represented by the LU genes, as well as calcium homeostasis and neuron differentiation represented by the DN genes (Fig. 5h; Additional file 2: Table S6). Moreover, we compared the nodes in our TRCs with cell type specific DEGs between no pathology and AD pathology obtained from single cell RNA sequencing of prefrontal cortex samples from 24 control individuals with no or little pathology (no-pathology) and 24 age-matched individuals with a spectrum of mild to severe β-amyloid and other pathologies (AD-pathology) [49]. The comparison revealed that Spi1 and Tnfrsf1 in the microglial TRC were shared with microglia specific DEGs identified from single cell RNA sequencing and Gadd45 in the OPC TRC was with OPC specific DEGs. Taken together, these data suggest that TRCs represent core regulatory programs underlying early and late alterations of cellular processes associated with prion disease.