Figure 1 shows the workflow of the integrative network analysis and functional validation experiments performed in this study. We first examined the relationships between GJA1 mRNA expression and clinical and pathological traits in 29 AD transcriptomic datasets. We also investigated the role of GJA1 in gene networks underlying AD. We then performed in vitro experiments to study the role of GJA1 in regulating AD gene networks and AD related phenotypes using primary astrocytes purified and cultured from wildtype and astrocyte specific Gja1−/− mice. Gja1’s target gene signatures identified from the RNA-seq data from the in vitro experiments were then projected onto the GJA1 centric networks to validate these networks’ structures.
GJA1 is a key regulator of an astrocyte specific gene subnetwork dysregulated in LOAD
Several studies have used co-expression network analysis to find modules of co-regulated genes in AD [36, 60, 61]. Our previous study developed a novel network approach capable of integrating clinical and neuropathological data with large-scale genetic and gene expression [98]. This network biology approach led to a novel multiscale network model of LOAD, which identified a number of coexpressed gene modules that were strongly associated with AD pathological traits or underwent dramatic disruption of high-order gene-gene interactions [98]. One such module, referred to as the khaki module in the original construction of this network, was of particular interest since it included APOE, the top AD risk factor gene. Moreover, the average interaction strength among its member genes in LOAD was reduced by 71% compared to that in normal control at a false discovery rate (FDR) < 2%, suggesting a huge loss of coordination among this group of genes in AD. The khaki module was enriched for the genes in Gamma-aminobutyrate (GABA) biosynthesis and metabolism (24 fold enrichment (FE), Fisher’s exact test (FET) p = 0.046) and harbored 12 (ALDOC, APOE, AQP4, ATP1A2, CSPG3, CST3, EDG1, EMX2, GJA1, PPAP2B, PRDX6 and SPARCL1) of 46 known astrocyte marker genes, a 15-fold enrichment over what would be expected by chance (FET p = 6.55E-9). The module was also enriched for the expression of the common variants identified as genome-wide significant by AD genome wide association studies (GWAS) (3-FE, FET p = 1.92E-11). Bayesian causal network analysis showed that GJA1 was the top driver of the module followed by FXYD1, STON2 and CST3 [98]. The key drivers of the corresponding causal network of the module were the nodes that had a large number of downstream nodes [90, 98]. These results indicate that GJA1 is a potential regulator of molecular networks in AD. In the next subsection, we will investigate the association between Gja1 mRNA expression and AD. All of the statistical significance levels reported were corrected for multiple testing unless otherwise specified.
GJA1 expression is associated with AD clinical and pathophysiological traits
To gain insight into the role of GJA1 in cognitive functions and AD pathogenesis, we first extensively investigated how GJA1 expression at the mRNA level was correlated with AD neuropathological traits in 29 gene expression datasets from three AD cohort studies of aging and dementia that included organ donation at death: the Mount Sinai/JJ Peters VA Medical Center Brain Bank (MSBB; Additional file 1: Table S1) [91], the Religious Orders Study and the Rush Memory and Aging Project (ROSMAP) [8, 9] and in the Harvard Brain Tissue Resource Center Alzheimer’s Disease study (HBTRC) [98]. We chose six different clinical and pathological criteria to evaluate the clinical relevance of GJA1 on AD pathology and cognitive functions: the MiniMental State Examination (MMSE) score [25, 32], the sum of NFT density estimates for all cortical regions examined (NTrSum), Mean Plaque density (PLQ_Mn) for the estimation of average plaque density, Braak stage score for quantitative assessment of neurofibrillary tangles [11], the Consortium to Establish a Registry for Alzheimer’s disease (CERAD: 1 for definite AD, 2 for probable AD, 3 possible AD, 4 for normal control) score for quantitative measure of neuritic plaques, and clinical dementia rating score (CDR ranging between 0 and 5 with 0 for normal control and 5 for severe dementia).
In the microarray data in the ROSMAP cohort, GJA1 expression was significantly correlated with CERAD score (r = − 0.15, p = 3.3E-3) and the MiniMental State Examination (MMSE) score (r = − 0.14, p = 6.4E-3). Similar results were observed in the ROSMAP RNA-seq dataset (Fig. 2a, and Additional file 1: Table S2), suggesting that the mRNA expression of GJA1 is associated with AD pathogenesis and dementia. The MSBB AD cohort includes microarray and RNA-seq data from a battery of distinct brain cortical regions and thus provides an excellent opportunity to investigate regional differences in the correlation between GJA1 expression and AD neuropathological traits [91]. Among the 19 brain cortex regions investigated in the MSBB AD microarray data, GJA1 expression in six cortex regions including BM10 (frontal pole), BM20 (inferior temporal gyrus), BM21 (middle temporal gyrus), BM32 (anterior cingulate), BM36 (parahippocampal gyrus) and BM46 (dorsolateral prefrontal cortex) was significantly correlated with at least three AD neuropathological traits (Fig. 2b, and Additional file 1: Table S2). Overall, GJA1 expression in these six cortex regions displayed a significant positive correlation with Braak stage score, PLQ_Mn, NTrSum and CDR. The MSBB AD RNA-seq data revealed a consistent pattern of correlation between GJA1 expression and AD clinic traits across the cortical regions studied (Fig. 2c-e; Additional file 1: Table S2). Notably, in BM10, BM36 and BM44 cortex regions, the microarray and RNA-seq data converged to show a consistent correlation between GJA1 expression and AD neuropathological traits (Additional file 1: Table S2). Thus, in the MSBB cohort, the association between GJA1 expression and AD neuropathological traits was cortex-specific. At the protein level, GJA1 in the brain cortex BM10 region was significantly correlated with CERAD (r = − 0.32, p = 1.26E-07), PLQ_Mn (r = 0.37, p = 4.08E-10) (Fig. 2f) and CDR (r = 0.35, p = 5.69E-9) (Additional file 1: Table S2). Also the total soluble amyloidβ (Aβ) levels had a significant positive correlation with GJA1 protein levels in the BM10 region (r = 0.18, p = 0.0036).
We showed previously in the HBTRC cohort that Gja1 had a significant correlation with Braak score (r = 0.61 and 0.52 for dorsolateral prefrontal cortex and cerebellum cortex regions, respectively) [98].
We further checked how GJA1 mRNA expression was correlated with 30 known AD risk factor genes. As shown in Fig. 2g, a majority of those AD genes were significantly correlated to GJA1 in the five RNA-seq datasets.
We also examined GJA1 differential expression between various subgroups of AD severity with respect to each individual AD neuropathological or functional cognitive trait using pairwise Student’s t-test. Consistent with the correlation analysis, Gja1 expression increased significantly as the disease deteriorated (Additional file 1: Table S3).
Because APOE is one of the major AD risk factors, and age and sex are critical clinic covariables in AD neuropathology, we investigated whether age, sex and APOE genotypes had any impact on the association of GJA1 expression with AD clinical and pathological traits. We stratified this cohort by age of death (AOD) (AOD > 85 versus AOD < 85), sex (female versus male), and APOE genotypes (E23, E34 and E33). As demonstrated in Additional file 1: Table S4 and Additional file 2: Figure S1, GJA1 expression is more significantly associated with clinical dementia rating (CDR) in the group with AOD > 85 than that with AOD < 85, in females than males, and in the group with APOE E33 than that with E34 or E23, across four brain regions in the MSBB cohort, suggesting that age, sex and APOE genotypes impact the association of GJA1 expression with clinical and pathological traits.
We further assessed if GJA1 mRNA expression was correlated with the variants of the known AD risk genes using the RNA-seq data from the brain region BM36 in the MSBB cohort. Among the 28 ADGWAS genes that had identifiable variants in the present study, 18 had at least one variant that possessed a significant correlation with GJA1 (Additional file 2: Figure S2). For example, the transcript/isoform ENST00000532146 of the risk factor CELF1 is significantly correlated with GJA1 expression in BM10 and BM44, but not BM22 and BM36 while ENST00000534614 and ENST00000539254 are correlated with GJA1 expression only in BM36 and BM22, respectively (Additional file 2: Figure S2). These results suggested that variation in transcripts’ abundance might be an important factor for determining the co-regulation between GJA1 and AD risk factors.
In summary, both correlation and differential expression analyses revealed that GJA1 was associated with amyloid and tau pathologies of AD as well as cognitive functions suggesting that GJA1 may play an important role in AD.
Transcriptomic changes caused by Gja1 deficiency in mouse astrocytes
To validate the role for GJA1 in orchestrating the astrocytic transcriptome, we purified and cultured primary astrocytes from wildtype and astrocyte specific Gja1−/− mice, and identified differentially expressed genes (DEGs) by RNA-seq in Gja1−/− vs wildtype primary astrocytic cultures in the absence or presence of wildtype primary cortical neurons. Each group had four replicates. All cultures were treated with 10 μM Aβ1–42 oligomer from div 10 through 14, when total RNAs were harvested. We identified 2891 upregulated (termed AST(up)) and 2605 downregulated (termed AST(dn)) DEGs upon the ablation of the Gja1 gene as compared to wildtype primary astrocytes (Fig. 3a). We identified 573 upregulated genes (termed AST + NEU(up)) and 1391 downregulated genes (termed AST + NEU(dn)) in Gja1−/− vs. wildtype primary astrocytes co-cultured with primary neurons (Fig. 3a). AST(up) shared 328 genes (11.3%) with AST + NEU(up) and 208 genes with AST + NEU(dn) while AST(dn) shared 672 and 60 genes with AST-NEU(dn) and AST-NEU(up), respectively (Fig. 3a). These DEG signatures were enriched for a variety of biological pathways including translational processes, immune response, cell-cell communication, extracellular matrix, microtubule cytoskeleton, synaptic transmission, lipid biosynthesis, steroid biosynthesis, cholesterol biosynthesis and cell-matrix adhesion (Fig. 3a-b, Additional file 1: Table S5).
AST(up) was significantly enriched for the genes in an Aβ network signature [16] (2.03 FE, FET p = 1.08 E-10), the AD genome-wide significant risk factor gene signature (ADGWAS) (3.47 FE, FET p = 6.58E-6) (Fig. 2c) while AST + NEU(dn) was enriched for the genes in the Aβ network signature (1.68 FE, FET p = 1.1E-3) (Fig. 3c and Additional file 1: Table S6). These results suggested that Gja1 is a critical regulator of the AD GWAS genes and may play an important role in Aβ metabolism.
We further intersected the GJA1 KO DEG signatures with the signatures from other inflammatory diseases to better understand of the immune response component in GJA1 regulated gene expression. We considered 2 well-established inflammatory gene signatures, an inflammatome signature of 2461 genes from eleven rodent inflammatory disease models [90] and the human macrophage and immune response enriched module of 2483 genes causally linked to obesity and diabetes [23]. These two inflammatory signatures share 758 genes [termed the core disease-related inflammatory gene (CDIG)]. About half of the CDIGs fall into AST(up) (accounting for about 12.6% of AST(up), 6.57 FE, FET p = 1.85E-208) and the intersection includes well known inflammatory markers (e.g., CD44, CD53, FCER1G, HCK, TYROBP and TREM2), inflammation complement component members (C1QA, C1QB, and C1QC), CXC chemokines (CXCL10, CXCL3 and CXCL6) and TNF Receptors (TNFRSF11B and TNFRSF13B). The CDIGs are also significantly enriched in AST + NEU(up) (accounting for about 15% of AST + NEU(up), 7.84 FE, FET p = 1.06E-49). These results suggest that inflammation is a critical component in GJA1-regulated gene expression.
We then investigated expression changes of astrocyte- and neuron-specific genes in various gene signatures using the recently identified brain cell type specific signatures [57]. Astrocyte- and neuron-specific marker gene signatures are enriched in the down- and up-regulated gene signatures in the neuron and Gja1−/− astrocyte co-cultures versus the neuron and wildtype astrocyte co-cultures, respectively (Additional file 1: Table S9A). The co-culture systems with and without Gja1−/− upregulated a significant portion of the neuron-specific marker genes when compared with the respective astrocyte alone models (with or without Gja1−/−) while down-regulating many astrocyte marker genes (Additional file 1: Table S9B). The gene lists from the abovementioned intersection analyses can be found in (Additional file 1: Tables S9C and S9D).
Gja1 deficiency induced transcriptomic changes highly overlap the GJA1 centric gene networks in AD
We further examined the molecular mechanisms of GJA1 in AD pathogenesis and cognitive function. We first examined if these Gja1−/− DEG signatures were enriched in the Khaki module where GJA1 resided based on our previous study [98]. As shown in Additional file 1: Table S7, AST + NEU(dn), AST(dn) and AST(up) were all significantly enriched in the khaki module with FET p = 3.74E-54 (9.0 FE), 3.64E-28 (4.4 FE), 1.28E-5 (2.1 FE), respectively. On the other hand, the overall signatures AST(all) and AST + NEU(all) also significantly overlapped the module with FET p = 3.03E-32 (3.2FE) and 1.9E-47 (6.7 FE), respectively .
To gain more insights into the signaling circuit of the GJA1 regulation in AD, we constructed the Bayesian causal network for the khaki module and projected the Gja1−/− DEG signatures onto this network in order to delineate the underlying causal relationship among the molecular constituents of this module. As shown in Fig. 4a and Additional file 1: Table S7, 13 of the 19 predicted key drivers except for GJA1 in the khaki module were regulated by GJA1 (4.94 FE, FET p = 4.70E-05). Strikingly, 11 out of the 13 GJA1 regulated key driver genes were within the GJA1’s downstream network neighborhood (Additional file 2: Figure S3A). To better gain insights into the functions of the Bayesian causal network in Fig. 4a, we built up a GIA1 signaling pathway map. As shown in Additional file 2: Figure S3B, GIA1 regulates a diversity of pathways such as regulation of gap junction activity, innate immune system, TGFβ signaling, DNA repair and lipid metabolism. These results suggested that Gja1 deficiency significantly impacted the khaki module and thus strongly validated our previous prediction of GJA1 as a driver of this important AD related subnetwork.
To more comprehensively determine GJA1’s target genes in AD, we further systematically identified the genes correlated with GJA1 in a number of human AD cohorts. Genes significantly correlated with GJA1 (Benjamini-Hochberg (BH)-corrected p-value < 0.05) were first identified in each of eight datasets including four from the MSBB RNA-seq cohort, one from the ROSMAP RNA-seq data and three from the HBTRC cohort. Note that only the data from the AD subjects were used for the correlation analysis. A majority voting method was utilized to calculate the frequency of each correlation across the eight datasets (Additional file 1: Table S8). We defined consensus GJA1-centered correlation signature (CGCCS) as a function of frequency threshold n, i.e., CGCCS(n) = {g | frequency(r(g, GJA1)) ≥ n}, where r(g, GJA1) represents a significant correlation between a gene g and GJA1 at FDR < 0.05, and n = 1, 2, …, 8. This process led to GJA1 centered correlation networks, also denoted as CGCCS(n).
Additional file 2: Figure S4 shows the enrichment of the GJA1 centered correlation networks for the previously identified DEGs in AST and AST + NEU. CGCCS (4) was most significantly enriched for the AST(all) and AST + NEU(all) DEG signatures with FET p = 1.0E-330 (3.2 FE) and 7.5E-311 (3.9 FE), respectively (Black lines, Additional file 2: Figure S4). CGCCS (6) is enriched for the DEG signature in AST with FET p = 6.2E-309 (3.8 FE). Figure 4b shows the network CGCCS (6) which includes 201 up-regulated (red nodes, termed GJA1_centered_AST(up)) and 307 down-regulated (blue nodes, termed GJA1_centered_AST(dn)) genes upon the knockout of Gja1 in astrocytes. About 30% of the genes in each of the two GJA1_centered_signatures belonged to the khaki module (FET p = 1.52E-90, 32.37 FE). This GJA1-centered correlation network included eight driver genes (AGT, BMPR1B, CST3, FXYD1, NTSR2, SLC15A2, SPON1, and STON2) in the khaki module [98]. Similar results were derived for the AST + NEU(up) and AST + NEU(down) signatures. Therefore, GJA1 impacts the expression of many key network drivers in the astrocytic subnetwork, suggesting its central role in this AD-related gene network. Furthermore, the analysis of the genes specific to the A1/A2 astrocytic activation in our DEG from Gja1−/− astrocytes [50] revealed that most of pan-, A1-, and A2-specific genes were generally upregulated (Additional file 2: Figure S5).
Inflammatory cytokines downregulate expression of Gja1 and the astrocytic subnetwork module
Previously it has been shown that LOAD-relevant inflammatory cytokines such as TNFα and IL-1β downregulated expression of Gja1 in astrocytes [19, 74]. Wildtype astrocytes were treated with TNFα, or IL-1β, or both for 7 days, confirming that Cx43 (Gja1 protein) was profoundly and synergistically reduced by both cytokines (Fig. 5a-b). Paradoxically, IL-1β significantly increased, but TNFα significantly decreased, Apoe protein levels (Fig. 5b). Following cytokine treatment, we tested Gja1, Apoe and eight key network driver genes found to be differentially expressed in Gja1−/− astrocytes by RNA-seq analysis, as a proxy to capture the network changes. qPCR analysis revealed that Gja1, Apoe and the other astrocytic subnetwork drivers were similarly downregulated by these cytokines (Fig. 5c).
Gja1 channel activity increases the expression of the astrocytic subnetwork module
Since these cytokines inhibit GJC and potentiate hemichannel activities [74], we asked whether inhibition of GJC and hemichannel activities regulates Gja1 and other astrocytic subnetwork drivers. Treatment of wildtype astrocytes with carbenoxolone (CBX, inhibitor of GJC and hemichannel [2, 95]) or lanthanum (La3+, inhibitor of hemichannel [2]) led to significant reduction of Gja1 and Apoe protein levels (Fig. 6a-b). Interestingly, the CBX treatment had broader effects on the reduction of the driver genes (Fig. 6c), while La3+ had generally milder and more selective effects, suggesting that GJC and hemichannel activities contribute to the regulation of distinct sets of the genes.
We next tested whether Gja1 channel activation can alter the astrocytic subnetwork. Gja1 hemichannel activity can be increased by quinine [81]. Treatment of wildtype astrocytes with quinine significantly upregulated Gja1 and Apoe protein and mRNA levels, along with transcriptional upregulation of a subset of the driver genes (Additional file 2: Figure S6). These data collectively supports Gja1, and specifically its channel activity, as an important regulator of astrocytic gene coexpression network including Apoe.
Increased neuronal survival following Aβ treatment when co-cultured with Gja1−/− astrocytes
We assessed the role for Gja1 in neuronal death and viability in a transwell co-culture system, in which cortical neurons were grown on the bottom of the plate, while astrocytes were cultured on the transwell insert. This allowed physical separation of the two cell types, while maintaining chemical continuity. Incubation of wildtype neurons and wildtype astrocytes with the indicated concentration of Aβ1–42 oligomers resulted in increased neuronal death [t (4) = − 2.941, p = 0.030 at 2 μM; t (4) = − 5.857, p = 0.004 at 20 μM] and decreased viability [t (6) = − 5.395, p = 0.002 at 2 μM; t (6) = 4.671, p = 0.003 at 20 μM] at day 7 (but not at day 4 (Fig. 7a)), in a dose dependent manner (Fig. 7b). Wildtype neurons cultured with Gja1−/− astrocytes were significantly resistant to death and loss of viability (Fig. 7b), consistent with a previous study showing that Gja1−/− deficient astrocytes conferred neuroprotection [69].
Reduced Apoe secretion and synthesis in Gja1−/− astrocytes
APOE is a well-known risk gene for LOAD and primarily produced by human and mouse astrocytes as well as microglia and neurons [52]. In addition, APOE was in the coexpression network module governed by GJA1. Therefore, we tested how Apoe levels were affected in the absence of Gja1 in mouse primary astrocytes. Quantitative PCR revealed significantly reduced Apoe expression in Gja1−/− astrocytes (t (6) = − 5.330, p = 0.002) (Fig. 8a). Consistently, we found lower Apoe protein levels in Gja1−/− astrocytes by immunoblot analysis of the lysates, with or without Aβ1–42 oligomer treatment (Fig. 8b). We analyzed Apoe in the conditioned medium from these astrocytes and found that secreted Apoe in the conditioned medium was consistently reduced in Gja1−/− astrocytes compared to wildtype, with or without Aβ1–42 oligomer treatment (Fig. 8c). Two way ANOVA revealed that there were significant main effects by genotype and treatment (F(1, 18)=121.618, p < 0.0005, F(2, 18)=9.364, p = 0.0016, respectively) but there was no interaction of genotype and treatment (F(2, 18) = 0.305, p = 0.741].
Neurons cocultured with Gja1−/− astrocytes produce higher levels of Aβ species
Since our RNA-seq analysis revealed that Gja1 regulated genes were highly enriched in Aβ network, we therefore analyzed Aβ production in neuron/astrocyte cocultures. Indeed, we observed higher levels of both Aβ40 and Aβ42 species in the conditioned medium from neurons co-cultured with Gja1−/− astrocytes at all time points examined [Aβ40; t(4) = 17.684, p < 0.0005 (div6), t(4) = 3.447, p = 0.026 (div13), t(4) = 11.112, p < 0.0005 (div16), t(4) = 11.112, p < 0.0005 (div20): Aβ42; t(4) = 8.422, p < 0.0005 (div6), t(4) = 2.047, p = 0.110 (div13), t(4) = 11.245, p < 0.0005 (div16), t(4) = 8.016, p < 0.001 (div20)] (Fig. 8d-e). Because astrocytes have been shown to take up Aβ, preferentially oligomers [67], we investigated the role for Gja1 in controlling Aβ clearance. When we tested whether wildtype and Gja1−/− astrocytes were proficient in taking up exogenously applied fluorescently labeled Aβ1–42 oligomers, we found that the number of Gja1−/− astrocytes with cell-associated Aβ1–42 oligomers was significantly reduced compared to wildtype astrocytes [t(20) = − 6.670, p < 0.0005] (Additional file 2: Figure S7).
Reduced neuronal activity in response to Aβ treatment when co-cultured with Gja1−/− astrocytes
In order to understand the impact of the network perturbation by loss of astrocytic Gja1−/− on neuronal functions, we co-cultured wildtype cortical neurons with either wildtype or Gja1−/− astrocytes on multielectrode arrays and compared spontaneous and stimulated neuronal activities before and after Aβ1–42 oligomer treatment (Fig. 9a). Although few significant differences were detected prior to Aβ treatment at 16 days in vitro (div), at 18 and/or 20 div (2 and 4 days after 10 μM Aβ1–42 treatment) most metrics of spontaneous activity (i.e. number of spikes, bursts, network bursts, and spikes per network burst) were significantly reduced in Gja1−/− astrocyte-cocultured neurons (Fig. 9b). The reduced spontaneous neuronal activity at div 20 suggested that co-culture with Gja1−/− astrocytes made neurons less responsive to Aβ-induced changes. While stimulated neuronal activity at div 14 was not affected by co-culture with Gja1−/− astrocytes (Fig. 10a), it was reduced at 20 div (number of spikes, spikes per burst, mean firing rate, and burst duration) (Fig. 10b-c).