Disrupted myelination network in the cingulate cortex of Parkinson's disease.

The cingulate cortex is part of the conserved limbic system, which is considered as a hub of emotional and cognitive control. Accumulating evidence suggested that involvement of the cingulate cortex is significant for cognitive impairment of Parkinson's disease (PD). However, mechanistic studies of the cingulate cortex in PD pathogenesis are limited. Here, transcriptomic and regulatory network analyses were conducted for the cingulate cortex in PD. Enrichment and clustering analyses showed that genes involved in regulation of membrane potential and glutamate receptor signalling pathway were upregulated. Importantly, myelin genes and the oligodendrocyte development pathways were markedly downregulated, indicating disrupted myelination in PD cingulate cortex. Cell-type-specific signatures revealed that myelinating oligodendrocytes were the major cell type damaged in the PD cingulate cortex. Furthermore, downregulation of myelination pathways in the cingulate cortex were shared and validated in another independent RNAseq cohort of dementia with Lewy bodies (DLB). In combination with ATACseq data, gene regulatory networks (GRNs) were further constructed for 32 transcription factors (TFs) and 466 target genes among differentially expressed genes (DEGs) using a tree-based machine learning algorithm. Several transcription factors, including Olig2, Sox8, Sox10, E2F1, and NKX6-2, were highlighted as key nodes in a sub-network, which control many overlapping downstream targets associated with myelin formation and gliogenesis. In addition, the authors have validated a subset of DEGs by qPCRs in two PD mouse models. Notably, seven of these genes,TOX3, NECAB2 NOS1, CAPN3, NR4A2, E2F1 and FOXP2, have been implicated previously in PD or neurodegeneration and are worthy of further studies as novel candidate genes. Together, our findings provide new insights into the role of remyelination as a promising new approach to treat PD after demyelination.

INTRODUCTION

Parkinson's disease (PD) is the second most common neurodegenerative disorder, clinically characterised by progressive deterioration of motor system due to the substantial loss of dopaminergic neuron in the substantia nigra pars compacta. Additional symptoms include postural instability and cognitive impairment. PD is a complex disease with many underlying causes. The majority of PD is sporadic. About 10% of PD patients are inherited with genetic causes or have a family history [1]. Genetic research has identified several key genes involved in PD pathogenesis, such as PINK1 [2], Parkin [3], DJ‐1 [4] and ATP13A2 [5] as autosomal recessive; meanwhile, SCNA [6, 7] and LRRK2 [8, 9] as autosomal dominant.

Functions of these PD genes have been intensively studied. For autosomal recessive PD, mutations in PINK1 and Parkin represent the most common causes [1]. PINK1 encodes the phosphatase and tensin homologue (PTEN)‐induced kinase 1, which is a serine‐threonine protein kinase with an N‐terminal mitochondria‐targeting sequence [2]. Subsequent studies have confirmed its kinase activity and revealed its subcellular distribution in the cytosol and mitochondria [10, 11]. Parkin is an E3 protein‐ubiquitin ligase [3]. Studies suggested that PINK1 and Parkin act in a common pathway in the clearance of damaged mitochondria via autophagy, a process known as mitophagy [12, 13, 14, 15]. PD‐related mutations in PINK1 and Parkin could disrupt the mitophagy process as shown in mammalian cells ([14]), which is consistent with the previous findings that mitochondrial abnormalities are potentially involved in the development of PD [16]. Furthermore, mutations in SCNA contributed to the autosomal dominant phenomenon for PD. SCNA encodes α‐synuclein, a soluble presynaptic protein [7]. Missense mutations, duplications and triplications of the entire SCNA gene have been reported and associated with PD, particularly for those early‐onset patients [17]. Genetic alterations of SCNA cause α‐synuclein aggregation and spread in the brain, which is a major constituent of Lewy bodies [7, 17]. Aggregates of misfolded α‐synuclein have been observed in sporadic PD patients, suggesting that α‐synuclein plays a central role in the pathogenesis of PD [18]. Transgenic mice over‐expressing either human wild type or mutant forms of α‐synuclein have been generated to investigate the pathological mechanisms of PD [19, 20, 21].

PD results from the death of dopaminergic neurons in the substantia nigra, striatum and other dopaminergic nuclei. Although PD began with defects of motor neurons, the majority of PD patients eventually showed different levels of dementia [22]. Recent efforts have revealed that cognitive impairment in PD involves dysfunction of the cingulate cortex. One study used cross‐sectional and longitudinal brain perfusion SPECT analyses to explore changes while developing dementia in PD. Their results showed that PD patients with dementia had significantly reduced perfusion in the right posterior cingulate, indicating a possible relationship between dementia in PD and cingulate cortex [23]. The dopaminergic innervation and the distribution of D1/D2 receptors in the cingulate cortex have been examined. PD patients with memory deficits showed more pronounced reduction in D2 receptor binding than PD patients with normal cognition [24]. However, how the cingulate cortex contributes to PD pathology is still largely unknown. At the transcriptional level, the disease‐associated changes in the cingulate cortex would provide valuable information to study the aetiology of PD. In this study, we analysed publicly available RNAseq data of the postmortem cingulate cortex from PD patients and healthy controls [25]. Our results indicated that glutamate receptor signalling and genes related to regulation of membrane potential were broadly upregulated, in agreement with previous studies showing involvement of abnormally increased glutamatergic transmission and membrane potential in postsynaptic neurons derived from PD patients [26, 27, 28]. Meanwhile, genes associated with myelination and oligodendrocyte development were dramatically downregulated, implying a dysfunctional myelination network in the progress of PD.

In the central nervous system (CNS), different cell types look wildly different, and carry out very different roles. We examined the behaviours of different CNS cell types in PD according to the established cell‐type‐specific gene expression signatures. We found that myelinating oligodendrocytes were the only cell type affected with statistical significance. More importantly, the downregulation of myelination pathways was also observed and shared in the cingulate cortex of DLB patients, indicating the similarity of underlying causes for cognitive deficits in Lewy body diseases.

Through integrated analysis with ATACseq data from human anterior cingulate cortex (ACC) [29], we inferred regulatory relationships between transcription factors and their possible targets for DEGs. We further divided the inferred GRN into two communities using the Louvain method, which contains a sub‐network important for myelin formation and oligodendrocyte development. Collectively, our data suggests that various pathways could be manipulated to improve treatments for PD. Moreover, a few DEGs were validated in two different PD mouse models: Pink1 knock‐out (Pink1−/−) and M83 transgenic mice—human α‐synuclein A53T transgenic mice driven by mouse prion protein (PrP) promotor (Tg M83). Further characterisation will be needed to elucidate their specific roles in PD.

MATERIALS AND METHODS

Animals

Pink1−/− mice generated in C57BL/6 genetic background was described as before [30]. Mice over‐expressing A53T‐SNCA (Tg M83) were obtained from the Jackson Laboratory (stock number 004479) [31]. Three male and one female Pink1−/− mice at the age of 6 months, and five male Tg M83 mice at the age of 8 months were used in qPCRs. Age‐matched C57BL/6 mice were used as qPCR control. Animals were housed and maintained under standard conditions with 12 h light/dark cycle, and free access to lab chow diet and water. All experiments and procedures involving animals were approved by the Animal Care and Use Committee of Wenzhou Medical University.

Cingulate cortex preparation and RNA extraction

Mice were anaesthetised with isoflurane and decapitated. The cingulate cortex was dissected out and prepared according to the Franklin and Paxinos mouse Atlas [32]. Briefly, the desired brain area was cut with DNase/RNase‐free scissors and blade, and then stored in RNAlater solution (Thermo Fisher Scientific). Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific).

RNAseq analysis

Differential expression analysis was performed with edgeR, limma, and DESeq2 [33, 34, 35]. Genes with less than 100 summarised counts for 12 samples were filtered out first. Genes are considered as differentially expressed when they had a Log2 (fold change) > Log21.5 or <−Log21.5, and p‐value < 0.05. The Venn diagram was generated using the VennDiagram package in R.

Functional enrichment analysis on the differentially expressed genes was conducted by using clusterProfiler package [36] to determine if the genes are enriched for specific terms. Enrichment tests for GO terms were calculated based on the hypergeometric distribution. Given a set of genes, the enrichment GO categories were returned after multiple hypothesis tests using the Benjamini–Hochberg method to correct the p‐values and control the false discovery rate (FDR).

ATACseq analysis

ATACseq data for human anterior cingulate cortex (ACC) was obtained from BOCA (Brain Open Chromatin Atlas) [29], which is a collection of maps of neuronal (NeuN+) and non‐neuronal (NeuN‐) chromatin accessibility across 14 distinct brain regions of five adult individuals.

ATACseq peaks were annotated by assigning the region to the nearest gene using ChIPseeker (tssRegion = c(−3000, 3000), overlap = “all”) [37]. Peak sequences for each annotated gene were extracted individually. Predicted TF binding sites within the peak sequences were determined using FIMO (q‐value < 0.1) [38].

Gene regulatory network construction and visualisation

Transcription factors (TFs) and their binding sites were determined based on the Animal Transcription Factor database (AnimalTFDB), HOCOMOCO and JASPAR databases. Regulatory relationships for DEGs were initially inferred using GENIE3, and then refined by TF binding sites information from ATACseq data. The GRNs were clustered into two major communities using Louvain modularity clustering method. The visualisation of the network was performed using Gephi (0.9.2).

qPCRs

RNA quality and quantity were assessed with a NanoDrop One (Thermo Fisher Scientific). For qPCRs, total RNAs were reverse transcribed with the HiScript II 1st Strand cDNA Synthesis Kit (Vazyme) according to the manufacturer's instructions. Real‐time PCR assays were performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme). GAPDH was used as a reference control. Primer sequences used in this study were listed in Table S1.

Statistics

All data were expressed as mean ± SEM. Data were analysed by a two‐tailed Student's t‐test between two groups (experimental and control). The criterion for statistical significance was p < 0.05.

RESULTS

DEGs identified with transcriptomic analysis

RNAseq data of the cingulate gyrus from eight healthy and eight PD patients were obtained from GEO with accession number GSE110716. The cingulate gyrus from PD patients were positive with an antibody against aggregated α‐synuclein [25]. The overall process of transcriptomic, regulatory network analysis and validation is illustrated in a flowchart (Figure 1a).

Venn diagrams of differentially expressed genes (DEGs) using DESeq2, edgeR, and limma (p < 0.05 and log2FC > log21.5 or <−log21.5) and functional enrichment analysis. (a) Flow chart of the RNAseq and regulatory network analysis and validations in this study. (b) Three programs were used for differential expression analysis of RNAseq data from the cingulate cortex of six PD patients versus six healthy controls. Intersection of the results, 889 DEGs, were used for downstream analysis. (c) Volcano plot of DEGs. NS: not significant. (d) Heatmap of top DEGs based on adjusted p values for eight PD patients versus eight healthy controls. (e) Heatmap of top DEGs based on adjusted p values for six PD patients versus six healthy controls. (f) Circular visualisation of GO enrichment and pathway analysis of all DEGs. Red dots represent upregulated DEGs, and blue dots represent downregulated DEGs. The outer circle shows annotated GO terms with associated DEGs. (g) Top 10 enriched GO terms in the category biological process for 255 up‐ and 402 down‐regulated DEGs, respectively. (h) Gene Concept Network plots (cnetplots) highlight the linkage between top enriched pathways such as demyelination and oligodendrocyte differentiation and associated DEGs. Colour intensity represents log2 fold changes of gene expression

Three different programs were used to determine the differentially expressed genes for the RNAseq data set: DESeq2, edgeR and limma, which are among the most cited transcriptomic analysis tools so far. As previous studies suggested, each package has its advantages and disadvantages, and no single method outperforms the other ones [39, 40]. Limma uses the linear models to analyse designed experiments and assess differential gene expression [34, 41]. In contrast, edgeR and DESeq2 use different normalisation approaches and dispersion estimates based on negative binomial distribution [33, 35]. As there are marked differences between the different algorithms, the three packages will return somewhat different sets of DEGs with considerable variability. The intersection of three sets represents a conservative estimate of DEGs, which will be optimal for downstream analysis.

Because genes with extremely low counts across all samples provided little evidence for differential expression rather than noise instead, genes with total counts less than 100 for all 16 samples were filtered out first. Any genes identified by all three programs were selected as DEGs with a Log2 (fold change) > Log21.5 or <−Log21.5 and p‐value < 0.05. In total, 501 DEGs were shared by the three programs. However, hierarchical clustering and heatmap using top‐ranked DEGs (sorted by adjusted p‐value) failed to separate the eight healthy controls from eight PD patients (Figure 1d). Obviously, two controls and two PD samples provided much noise and outliers, which are likely due to common postmortem RNA degradation. After removing these four samples, random errors decreased, therefore, we were able to increase the statistical power. As a result, 889 DEGs were successfully identified by all three RNAseq packages, including 351 up‐ and 538 down‐regulated genes (Figure 1b, 1c and Table S2). PD patients and healthy controls formed separated clusters in the heatmap based on the top‐ranked DEGs (sorted by adjusted p‐values; Figure 1e). Subsequent analyses were performed using these 889 DEGs.

Functional enrichment analysis of DEGs

The pathology of Parkinson's disease would induce dramatic changes of gene expression in specific cellular and biological processes. First, gene ontology analysis was performed to find over‐represented functional categories for all DEGs. The most enriched biological processes include myelination, axon ensheathment, neuropeptide signalling pathway, response to stress, and regulation of membrane potential. Notably, expression of DEGs present in the following GO categories “GO:0042552: myelination”, “GO:0007272: ensheathment of neurons”, and “GO:0008366: axon ensheathment” were downregulated (Figure 1f).

Next, we ran enrichment analysis separately for up‐ and down‐regulated genes. Upregulated DEGs are related to modulation of membrane potential and synaptic transmission (Figure 1g). The most significant GO cellular component terms identified for upregulated genes are “synaptic membrane”, “transporter complex”, and “ion channel complex” (See Figure S1a). Major GO molecular function terms for upregulated genes include “channel activity”, “transporter activity”, and “neurotransmitter receptor activity” (See Figure S1b). Whereas, the downregulated DEGs are significantly enriched in following biological processes, such as “ensheathment of neurons”, “myelination” and “oligodendrocyte differentiation” (Figure 2a), indicating glial defects and demyelination in PD cingulate cortex. Furthermore, these pathways and their associated DEGs were visualised in a cnetplot (Figure 1h). Cellular component analysis showed that these downregulated DEGs are expressed in “apical plasma membrane”, “cell‐cell junction”, “main axon”, and “myelin sheath” (See Figure S1a). GO molecular function analysis revealed that some of the downregulated DEGs are related to “actin binding” (See Figure S1b). Previous studies have demonstrated that actin cytoskeleton plays a central role in myelination in both central and peripheral nervous systems [42, 43]. Together, functional enrichment analysis showed consistent results, suggesting demyelination could be a hallmark of the cingulate cortex in PD patients.

FIGURE 2

Expression changes of cell‐type‐specific signature genes in the Parkinson’s disease (PD) cingulate cortex and shared downregulated myelination pathways in the DLB ACC. (a) Box plots for the log10 ratios of gene expression levels (normalised counts) of cell‐type‐specific signature genes for PD patients and healthy controls were made in eight different CNS cell types. (b) Top 10 enriched GO terms in the category biological process for differentially expressed genes (DEGs) in DLB ACC. (c) 16 genes are shared between 538 downregulated DEGs of the PD cingulate cortex and 31 DEGs of the DLB ACC. (d) Top 10 enriched GO terms in the category biological process for the shared 16 genes

Cell‐type‐specific changes in the PD cingulate cortex

The bulk RNAseq dataset analysed in this study were from whole tissue samples of the cingulate cortex consisting of heterogeneous cell types. Transcriptomic signature changes in bulk brain tissues are significantly confounded by underlying differences in cell‐type composition. Therefore, it is crucial to dissect cell heterogeneity to unveil regulatory changes for disease aetiology.

To discriminate distinct cell types in the CNS, RNAseq studies have been used to generate a transcriptome database for 8 cell types from mouse cerebral cortex, including neurons, astrocytes, oligodendrocyte precursor cells (OPCs), newly formed oligodendrocytes (NFO), myelinating oligodendrocytes (MO), microglia, endothelial cells, and pericytes [44]. As a result, RNAseq analysis revealed distinct gene expression signatures for each cell type. We used a panel of 50 genes with the top 40 cell‐specific markers and top 10 transcription factors that have the highest mutual information for different cell clusters [44].

By using the reference panel, we examined overall changes of 8 CNS cell types. Table S3 summarised the gene lists for each cell type and their normalised expression counts. The results clearly showed that myelinating oligodendrocyte genes were significantly downregulated (p = 0.00073), while newly formed oligodendrocyte genes were also downregulated with no statistical significance (p = 0.18). Expression of both neurons and microglia genes appeared to increase in PD, however, not statistically significant (p = 0.43 and p = 0.16, respectively) (Figure 2a and Table 1). In summary, at the resolution of individual cell types, we identified that myelinating oligodendrocytes are a major cell type damaged in the PD cingulate cortex.

TABLE 1

Means of log10‐transformed normalised counts of signature genes used for each cell type in PD and healthy controls and the p‐values

	Astrocytes	Neurons	OPCs	NFO	MO	Microglia	Endothelial	Pericytes
CT	2.619737	2.217234	2.416256	2.645571	2.75292	1.746067	2.133549	2.097207
PD	2.614914	2.289216	2.427939	2.560417	2.512567	1.842782	2.073687	2.041602
p‐value	0.94	0.43	0.83	0.18	0.00073	0.16	0.26	0.38

Abbreviations: MO, myelinating oligodendrocytes; NFO, newly formed oligodendrocytes; OPCs, oligodendrocyte precursor cells; PD, Parkinson’s disease.

Validation of demyelination pathways in ACC of dementia with Lewy bodies (DLB)

All patients in the PD RNAseq dataset analysed above were positive for Lewy bodies in their cingulate cortex. Spread and aggregation of the α‐synuclein form structures called Lewy bodies (LB) in neuronal cell bodies and neurites in neuronal processes. Previous studies have shown a positive correlation between Lewy body densities in cortical regions including the cingulate cortex and cognitive impairment in PD [45, 46]. Parkinson's disease with dementia (PDD) and dementia with Lewy bodies (DLB), known as the Lewy body diseases, distinguished arbitrarily based on the timing of cognitive decline relative to motor symptoms, share important clinical features and common underlying molecular pathogenesis [47, 48]. Since dementia is common in PD patients, approximately 50%–80% PD patients will develop cognitive decline over the years [49], PD, PDD and DLB are viewed as a different spectrum of the same disease [47].

An RNAseq study of ACC reported for eight DLB patients and ten age‐matched controls, 490 DEGs were retrieved, 123 upregulated and 367 downregulated [48]. GO analysis indicated that the DEGs are mainly enriched in myelination, axon ensheathment and gliogenesis pathways (Figure 2b), which are consistent with the PD RNAseq analysis we conducted above. In total, 31 DEGs were extracted from these enriched GO terms. Sixteen overlapping genes were identified when intersected with downregulated 538 DEGs from our PD RNAseq analysis (Figure 2c). Top scored GO terms of the 16 overlapping genes are “gliogenesis”, “myelination”, and “ensheathment of neurons” (Figure 2d).

In conclusion, the downregulation of myelin genes and pathways in the cingulate cortex for both PD with Lewy bodies positive in the cingulate cortex and DLB could suggest a strong association between demyelination and cognitive disorders in these Lewy body diseased patients.

Gene regulatory network construction

Transcriptional regulatory changes are prominent features and mechanisms of human diseases. Neurodegenerative diseases such as AD and PD involved expression changes of thousands of genes in the brain. Convergent evidence has shown that dysregulation of TFs and their target genes contributed to disease risks and aetiology. Therefore, it is important to understand transcriptional regulation through network analysis.

RNAseq alone cannot determine direct interactions and often produced false positive interactions, the use of ATACseq data not only can validate direct TF‐target relationships based on TF binding sites information, but also help identify potential novel regulatory relationships. However, there is no matched ATACseq data for PD cingulate cortex available, instead we used ATACseq data from human anterior cingulate cortex (ACC) for assessment of chromatin accessibility and validation of TF binding sites in the target genes. The ATACseq data for human anterior cingulate cortex (ACC) was obtained from BOCA (Brain Open Chromatin Atlas) [29], which is a collection of maps of neuronal (NeuN+) and non‐neuronal (NeuN‐) chromatin accessibility across 14 distinct brain regions of five adult individuals. During the past several decades, cell counting studies have found that adult human brain contains about equal numbers of neuronal and non‐neuronal cells, and non‐neuronal cells are primarily glia [50, 51]. With three‐dimensional re‐construction of 50,000 cells from the temporal lobe of the human cerebral cortex, a new study found that glia outnumbered neurons at a ratio of 2:1 and oligodendrocytes were the most common non‐neuronal cell type [52].

For neuronal cells of ACC, ATACseq peaks were significantly enriched in enhancers, including distal intergenic (27.1%) and other introns (31.14%), compared to the proximal promoter (≤1 kb) (19.33%). In contrast, for non‐neuronal cells of ACC, the percentage of ATACseq peaks in proximal promoter (≤1 kb) were much higher, 32.51%; meantime, the percentages of ATACseq peaks located in enhancer regions were lower: 20.71% for distal intergenic and 25.21% for other introns. The difference in peak distribution might be due to different transcriptional regulations for neuronal and non‐neuronal cells/glia (Figure 3a). We further examined the distribution of ATACseq peaks for the downregulated myelin genes. Overall, there are stronger and more peaks in non‐neuronal than neuronal cells for these myelin‐related DEGs. For example, MYRF (Myelin Regulatory Factor), CD9 and NKX6‐2 have distinct non‐neuronal peaks present in their promoters and introns (Figure 3b). E2F1, a transcription factor regulating cell cycle progression, plays important roles for neurogenesis and myelination, which has similar peaks in both neuronal and non‐neuronal cells [53] (Figure 3b).

FIGURE 3

Chromatin accessibility analysis using ATACseq data in human ACC. (a) Proportions of ATACseq peak regions identified in neuronal and non‐neuronal cells that represent various genomic annotations. (b) Three representative myelin and oligodendrocyte genes, MYRF, CD9, and NKX6‐2, show distinct ATACseq peak distribution in non‐neuronal cells; E2F1 shows similar ATACseq peaks in neuronal and non‐neuronal cells

Gene regulatory networks (GRNs) capture the regulatory relationships between transcription factors (TFs) and their target genes. Many network inference methods based on gene expression data have been developed. GENIE3—a Random Forest‐based algorithm has consistently been one of the best performers [54, 55]. We applied GENIE3 to determine transcriptional regulatory relationships for the cingulate cortex of PD. Initially, 45 putative transcription factors (TFs) within DEG list were identified according to the Animal Transcription Factor DataBase (AnimalTFDB) [56]. The input of GENIE3 is an expression matrix of DEGs and a list of the 45 TFs, which were regarded as the candidate regulators. Each edge has a weight assigned by GENIE3, as a measurement of the confidence of the prediction. After regulator → target weighted edges produced and ranked by GENIE3, only edges with the top 25% of weights were retained for further analysis. Thus, we built preliminary TF‐target links based on gene expression data.

Out of the 45 putative TFs, 32 have binding sites in HOCOMOCO and JASPAR databases, two widely used TF binding profiles [57, 58]. Therefore, we were able to narrow down the number of TFs to 32. ATACseq peaks were annotated by assigning the region to the nearest gene using ChIPseeker [37]. Peak sequences for each annotated gene were extracted individually. Predicted TF binding sites within the peak sequences were determined using FIMO (q‐value < 0.1) [38].

The 32 TFs and their 466 target genes were inferred from the ATACseq data. Finally, we concluded GRNs with 1128 links for DEGs by intersecting the TF‐target links from ATACseq data and the TF‐target links from RNAseq data (See Table S4).

Network clustering and functional enrichment analysis

Community detection is necessary to understand the structure of large and complex networks. Therefore, after constructing GRNs, we ran a clustering analysis to find distinct clusters with high modularity within DEGs. We used the Louvain method, one of the most popular algorithms for uncovering hierarchies of communities [59]. The modularity measures the density of connections within clusters compared to the density of connections between clusters. Therefore, we were able to accurately identify sub‐ or sub‐sub‐networks within DEGs [59].

Python–Louvain algorithm produced two distinct clusters: one smaller, tightly grouped cluster—cluster 1 and another, large, dispersed cluster—cluster 2 (See Table S5). After clustering analysis, we used Gephi to visualise the network with the Force Atlas2 layout algorithm (Figure 4a). Interestingly, DEGs in cluster 1 are mainly upregulated, and only 5 out of 192 are downregulated; in contrast, cluster 2 contains 289 downregulated and 17 upregulated DEGs (Figure 4b).

FIGURE 4

Integrative analysis of RNAseq and ATACseq and modularity clustering of networks. (a) Network topology visualised through Gephi software. (b) Upset graph shows the distribution of up‐ and down‐regulated DEGs in two clusters. (c) Sub‐network for myelin and oligodendrocyte development. (d) Enrichment analysis for the two clusters

To better understand regulatory mechanisms for oligodendrocyte development and myelination in the cingulate cortex, we built a sub‐network for 18 DEGs involved in myelination and gliogenesis: 5 TFs and 13 target genes (Figure 4c). It is shown that there are reciprocal interactions between these 5 TFs and a target gene can be regulated by multiple TFs. The size of the circle represents the number of targets. Clearly, Olig2, Sox8, Sox10, E2F1 and NKX6‐2 are the most important regulators for myelination and gliogenesis. More importantly, functional enrichment analysis for these two clusters showed essentially similar results with up‐ and down‐regulated DEGs (Figure 4d). The top enriched GO terms for cluster 1 include regulation of membrane potential and glutamate receptor signalling pathway. DEGs in cluster 2 are significantly enriched in oligodendrocyte differentiation and myelination. In summary, through network and clustering analysis, at transcriptional level we discovered how neurons and glia in the cingulate cortex responded to PD. The results would advance our understanding of dysfunctional myelination pathway and abnormal synaptic transmission during PD pathogenesis.

Validation of DEGs with two PD mouse models

The use of animal models is critical to elucidate the disease mechanisms for biomedical research.

Current PD animal models have clear limitations to model PD patients, because none of the PD models could completely recapitulate the key neuropathological and clinical features of PD. Although mouse models of PD did not show degeneration of nigrostriatal DA system, we still have gained great insights into underlying pathogenic mechanisms of PD from collectively studying these existing models that do mimic many aspects of PD. Thus, it is important to validate novel PD‐related DEGs in various animal models and dissect out their functions.

We next sought to selectively validate DEGs in two PD mouse models: Pink1 KO (Pink1−/−) and M83 transgenic mice—human α‐synuclein A53T transgenic mice driven by mouse prion protein (PrP) promotor (Tg M83). We selected 30 PD or neurodegeneration‐related DEGs for qPCR validation in these two PD mouse models. Two of them were validated in both models, TOX3 and NECAB2, which were upregulated in PD cingulate cortex. TOX3 is a member of the TOX High Mobility Group Box Family, in which SNPs are implicated in breast cancer, restless legs syndrome, and PD [60]. NECAB2 (N‐Terminal EF‐Hand Calcium Binding Protein 2) is essential for calcium homoeostasis. In iPSC‐derived midbrain dopaminergic neurons from GBA1‐associated PD patients, the expression of NECAB2 is significantly increased and calcium homoeostasis is dysregulated [61]. Several other DEGs, NOS1, CAPN3 and NR4A2, were validated in Pink1−/− mice; E2F1 and FOXP2 were confirmed in Tg M83 mice (Figure 5 and Table 2). According to RNAseq results and previous literature, additional 29 DEGs that play important roles in the CNS, which have not been reported associated with PD, were also selected for validation. PRLR, COX7A1 and KCNG1 were confirmed in Tg M83 mice, AMIGO2, validated in Pink1−/− mice (Figure 6 and Table 3). However, DEGs involved in myelination and gliogenesis pathways were largely not validated in our QPCRs. The key functions of selected 59 DEGs and their validation status were summarised in Tables 1 and 2. Our results reminded us of the limitation of PD mouse models. Since a variety of PD animal models were generated to study different aspects of the disease, it is not surprising to find different sets of DEGs were validated in the two PD mouse models, respectively. Nevertheless, our validations provide a very useful resource for functional dissection of the molecular mechanism in the cingulate cortex responding to PD using experimental animals.

FIGURE 5

QPCR results for 30 Parkinson’s disease (PD)‐related differentially expressed genes (DEGs), and myelination and gliogenesis related DEGs are underlined. (a) DEGs were confirmed in Pink−/− mice. (b) DEGs were confirmed in Tg M83 mice. (c) QPCR results for myelination and gliogenesis‐related DEGs in Pink−/− mice. (d) QPCR results for myelination and gliogenesis‐related DEGs in Tg M83 mice. (e) QPCR results for the remaining DEGs in Pink−/− mice. (f) QPCR results for the remaining DEGs in Tg M83 mice

TABLE 2

qPCR results for selected differentially expressed genes (DEGs), which are functionally related to PD and neurodegenerative diseases. The underlined genes are associated with myelination and gliogenesis

Genes	Log2FC	qPCRs in		PD‐related functions	References
	PD/ctrl	PINK1 −/−	SNCA tg
TOX3	0.9056	↑	↑	Its variants are associated with PD.	[60]
NECAB2	0.7527	↑	↑	It might be involved in calcium homoeostasis in PD.	[61]
NOS1	0.6805	↑	NS	Its polymorphisms are associated with PD. NOS inhibitors relieved the symptoms of PD rats.	[96, 97]
RGS14	0.6569	NS	NS	Diminished expression and altered cellular distribution of RGS14 were observed in dopaminergic neurons from PD substantia nigra.	[98]
GPC4	0.6487	↓	NS	Reducing GPC4 expression increased the production of dopaminergic neurons and improved motor behaviour in a rat model of PD.	[99]
E2F1	−0.6854	NS	↓	E2F1 was activated in mature dopaminergic neurons in PD patients and mediated neuronal death.	[53, 100]
				E2F1 is a key regulator of proliferation/differentiation transition for oligodendrocyte progenitor cells.
LRP2	−1.5611	NS	NS	LRP2 is required for oligodendrocyte and astrocyte development.	[101, 102]
				Expression of LRP2 in PD nigra and frontal cortex were up‐regulated.
UGT8	−0.7844	NS	NS	Plays protection roles for ageing brain and implicated in regulatory networks of PD	[103]
OPALIN	−1.5076	NS	NS	In LRRK2 KO mice, both SYT2 and OPALIN expression in striatal tissue were significantly up‐regulated.	[104]
MOBP	−1.6722	↑	↑	Associated with PD by GWAS; in addition, MOBP is present in Lewy bodies in PD	[105, 106]
MBP	−1.7839	NS	↑	In α‐syn tg mouse models, the striatum exhibited an age and disease‐dependent loss of MBP, indicating of the loss of myelin.	[107]
FRMD4B	−0.8897	↑	NS	As a scaffold protein, it was identified as a prominent hub gene of myelination network in AD.	[108]
MOG	−1.3195	↑	↑	Expressed on the oligodendrocyte cell surface and the outermost surface of myelin sheaths, in MPTP‐induced mouse model of PD, MOG immunisation prior to MPTP administration significantly diminished the microglia activation and inflammatory reaction.	[109]
NXPH3	−0.6973	NS	↑	A synapse‐related peptide, supported the survival of the grafted dopaminergic neurons in vitro and in vivo.	[110]
DLL1	−1.0789	NS	↑	LRKK2 could modulate the notch signalling by regulating DLL1 trafficking, which in turn may alter the neurogenesis of dopaminergic neurons in PD.	[111]
ERBB3	−1.0582	NS	NS	ERBB1‐4 and EGF have extensive cross talk with neuregulins. Defects of their interactions have been implicated in the pathology of PD.	[112]
GPR37	−0.8127	↑	↑	A substrate of parkin	[113]
ABCA8A	−0.8180	↑	↑	Human ABCA8 has two orthologs in mouse genome: ABCA8A and ABCA8B.	[114]
	ABCA8			Overexpression of ABCA8 in oligodendrocytes promoted ectopic expression of α‐syn.
ABCA8B	−0.8180	NS	↑	See above	See above
	ABCA8
DAAM2	−0.8302	NS	↑	Methamphetamine‐induced release of dopamine and down‐regulation of DAAM2 on the injection side in hemi‐parkinsonian rats.	[115]
ONECUT2	−2.2924	NS	↑	A homoeobox transcription factor, it is involved in midbrain dopaminergic neuron development.	[116]
NR4A2	−0.9859	↓	NS	It is required for the development and maintenance of dopaminergic neurons; and it is associated with pathogenesis of PD.	[117, 118]
GJB1	−1.2142	↑	↑	GJB1 had a preferential uptake of α‐syn oligomeric assemblies in neurons and oligodendrocytes. Also, α‐syn accumulation was positively correlated with upregulation of GJB1.	[119]
HAPLN2	−1.2718	↑	↑	Hapln2 might promote the aggregation of α‐synuclein in PD.	[120]
SYT2	−1.3775	NS	NS	SYT2 was significantly lost in the prefrontal cortex of PD patients with dementia according to proteomics.	[121]
BCAS1	−1.4747	NS	↑	The downregulation of BCAS1 was observed for human PD postmortem substantia nigra and the subthalamic nucleus isolated from the MPTP mouse model of PD.	[122, 123]
CAPN3	−0.8295	↓	↑	There exists a positive feedback loop between calpain activation and Aβ aggregates in AD: Aβ enhances calpain activation, which in turn produces more Aβ and therefore more calpain activation.	[124]
EPB41	−0.7289	NS	NS	It is normally associated with cytoskeleton, recently identified as a novel candidate gene in neurodegenerative diseases by integrated meta‐analysis.	[125]
FOXP2	−0.8102	NS	↓	FOXP2‐driven network is implicated in neurodegeneration.	[126]
SLCO1A2	−1.1985	NS	↑	Solute carrier organic anion transporter family member 1A2, it is identified as a new susceptibility locus for parkinsonian neurodegenerative tauopathy.	[127]

FIGURE 6

QPCRs results for 29 other differentially expressed genes (DEGs), and myelination and gliogenesis related DEGs are underlined. (a) DEGs were confirmed in Pink−/− mice. (b) QPCR results for myelination and gliogenesis‐related DEGs in Pink−/− mice. (c) DEGs were confirmed in Tg M83 mice. (d) QPCR results for myelination and gliogenesis‐related DEGs in Tg M83 mice. (e) QPCR results for the remaining DEGs in Pink−/− mice. (f) QPCR results for the remaining DEGs in Tg M83 mice

TABLE 3

qPCR results for other selected differentially expressed genes (DEGs), which have unknown functions in neurodegeneration. The underlined genes are associated with myelination and gliogenesis

Genes	Log2FC	qPCRs in		Nervous‐system related functions	References
	PD/ctrl	PINK1 −/−	SNCA tg
PRLR	1.6991	NS	↑	Prolactin receptor, it is important for neuroprotection, neurogenesis and sex‐specific pain regulation.	[128, 129]
KCNG1	1.2455	NS	↑	It is involved in slower, sustained or rhythmic activation in the cortex.	[130]
AMIGO2	1.0124	↑	NS	It is a novel membrane anchor of PDK1 and activate PI3K/Akt signalling.	[131]
RGS8	0.9767	↓	NS	RGS8 is highly expressed in the brain and negatively regulates the MCH1R signalling.	[132]
GRID2	0.6647	NS	↓	It encodes the glutamate receptor subunit delta‐2 and its deletions cause cerebellar ataxia and eye movement abnormalities.	[133]
CHRDL1	0.6647	NS	↓	Its mutations lead to X‐linked megalocornea. Astrocyte‐secreted Chrdl1 promotes GluA2‐dependent synapse maturation and limits synaptic plasticity.	[134, 135]
MAL	−0.8081	NS	NS	MAL is a raft‐associated membrane protein involved in myelination.	[136]
TMEM63A	−0.8820	↑	↑	A mechanosensitive ion channel, highly expressed in oligodendrocytes. Its heterozygous variants have been implicated in myelination deficit.	[137, 138]
ANLN	−0.9067	NS	↑	It is a cytoskeletal adaptor protein, highly expressed in myelinating oligodendrocytes and is required for proper myelin septin assembly.	[139]
NDRG1	−0.9672	NS	↑	Its mutations caused hereditary motor and sensory neuropathy. Its expression in the cytoplasm of Schwann cells is essential for the maintenance of myelin sheaths in peripheral nerves.	[140]
CLDN11	−0.9967	NS	↑	It is a major component of CNS myelin and have essential structure functions in maintaining normal myelin.	[141]
PLP1	−1.0052	NS	↑	It is the most abundant myelin protein in the CNS. Axonal swellings and degeneration were observed in PLP1 KO mice.	[142]
MYRF	−1.2770	↑	↑	It is a transcription regulator that specifically activates and maintains the expression of myelin genes in oligodendrocytes.	[143]
ERMN	−1.2970	↑	↑	It is necessary for myelin sheaths formation and maintenance. Knock‐out of ERMN accelerated the demyelination process.	[144]
SH3TC2	−1.4747	↑	↑	It is required for myelination by interacting with Neuregulin‐1/ErbB pathway.	[145]
NKX6‐2	−1.7785	↑	↑	A transcription repressor, it regulates oligodendrocyte gene expression and myelination. Its mutations are associated with spastic ataxia and hypomyelination.	[146]
RHOU	−0.6924	↑	NS	An atypical Rho GTPase, implicated in a variety of cellular processes: Cell shape, adhesion, mobility and survival	[147]
PIP4K2A	−0.7116	↑	NS	PIP4K2A can act as a negative regulator of PI3K, which is essential for neuronal protection.	[148]
NDE1	−0.7850	NS	NS	Its mutations are associated with intellectual deficiency and schizophrenia.	[149]
SCN1B	−0.8552	NS	↑	Sodium channel β1, its mutation variants have been linked to childhood epilepsy.	[150]
NEK3	−0.9112	↑	↑	Its mutations caused abnormal neuronal morphology and polarity, indicating NEK3 could involve in axonal projection and degeneration.	[151]
FAM107B	−0.9420	NS	↑	In the adult mice, its expression is restricted to the dentate gyrus, but function is unknown.	[152]
COX7A1	−0.9796	NS	↓	It encodes cytochrome c oxidase polypeptide 7A1, related to energy metabolism.
PDE8A	−0.7402	NS	↑	Implicated in inflammatory cell activation, PDE8A was identified as the predictor of suicidal behaviour.	[153]
KCNH8	−1.0997	↑	NS	A potassium voltage‐gated channel, its function is largely unknown.	[154]
SH3RF2	−1.2437	NS	NS	Heterozygous mutant mice of SH3RF2 displayed deficits of hippocampal neurons and autistic‐like behaviours.	[155]
TMEM235	−1.8490	NS	NS	It might be related to tight junction proteins—claudins or the voltage‐dependent calcium channels.	[156]
ADAMTS14	−1.8903	↑	NS	Related to inflammation and immune response	[157]
ADAMTS18	−2.0856	NS	↑	Related to inflammation and immune response	[158]

DISCUSSION

In addition to motor deficits, non‐motor symptoms occur in up to 90% of PD patients, including cognitive impairment, anxiety, depression, sensory dysfunction, sleep disorders, and vision problems. There is increasing evidence that non‐motor symptoms of PD are not just due to the loss of DA neurons in the basal ganglia [62]. Research on brain circuits involving in other brain nuclei will be especially important to develop treatments for PD non‐motor features.

The cingulate cortex is an extensive area of limbic system, important for emotion, memory, and action‐outcome learning [63, 64]. Fluorodeoxyglucose (FDG) positron emission tomography (PET) in dementia imaging showed hypometabolism in the inferior parietal lobe and posterior cingulate/precuneus can be a predictor of cognitive decline from mild cognitive impairment to AD dementia [65]. Particularly, a recent MRI study in 159 PD patients showed that cognitive impairment and excessive daytime sleepiness were associated with loss of integrity and atrophy in the anterior and posterior cingulate networks [66]. Hence, alterations in the cingulate cortex might be crucial in the development of PD dementia.

In this study, we conducted transcriptomic and network analyses for the cingulate cortex of PD patients. Genes related to regulating membrane potential and glutamate receptor signalling pathway were upregulated. It has been determined that the glutamate concentrations in plasma and serum are higher in PD patients than in healthy subjects [27, 67]. Alterations of glutamate receptors, such as NMDARs and AMPARs, have been observed in PD patients and in experimental models as well. Preclinical and clinical studies suggest that compounds targeting glutamate receptors could ameliorate PD motor symptoms [26, 27]. In line with these previous studies, for the first time we uncovered broadly increased expression of glutamate receptors and synaptic transmission pathway genes in PD cingulate cortex. Furthermore, key oligodendrocytes and myelin genes, including transcription factors that control oligodendrocyte differentiation and myelin formation, were significantly downregulated, which is likely an indication of demyelination.

Remarkably, based on cell type gene signature panel, we identified myelinating oligodendrocytes as the only cell type significantly affected in the PD cingulate cortex. The dysregulation of myelination pathways was largely shared and validated in the ACC of DLB patients. Identification of alterations at the resolution of individual cell types will not only help us elucidate mechanisms of PD, but also develop tools that allow cell‐type‐specific genetic manipulations.

In the CNS, oligodendrocytes wrap around axons to create myelin sheath. Schwann cells produce myelin in the peripheral nervous system. Myelin is vital to a healthy nervous system, mediating rapid action potential transmission and providing trophic support of long axons. Therefore, myelination is essential for formation and synchrony of neural networks. Reductions in myelination could disrupt efficient nerve signal transduction and contribute to axon degeneration. Myelin levels decrease during normal ageing. Myelin deficits have been consistently observed in schizophrenia, depression and bipolar disorder [68, 69, 70, 71]. Besides multiple sclerosis (MS), the most common demyelinating disease, oligodendrocytes/myelination have attracted considerable attention for other neurodegenerative diseases. Pioneering studies have shown association between myelin breakdown and progression of AD and Huntington's disease [72, 73]. Additionally, using histological and immunostaining methods, demyelination and oligodendrocyte loss have been observed in the frontal lobe white matter and cortex grey matter in AD patients [74, 75]. Later, transcriptomic analysis revealed downregulation of myelination networks in AD and progressive supranuclear palsy [76]. However, studies on myelin changes in PD are scarce. A cross‐sectional neuroimaging study has implicated alterations of myelin content in PD brains [77]. Rather than being passive supporting cells as long thought, myelin/oligodendrocytes play an active role in brain function and information processing. Studies demonstrated generation of new oligodendrocytes and formation of new myelin are required for motor skill learning, consolidation, and retrieval of remote fear memory [78, 79, 80]. Thus, demyelination might contribute to memory loss and cognitive impairment in both AD and PD. Compared to AD, demyelination and myelin damage have been overlooked in PD; our results would add new and compelling evidence of demyelination in the progress of PD.

As an indispensable and essential component of CNS, researchers have shown remyelination can be a promising and effective therapeutic strategy [81]. In mouse models of multiple sclerosis, two drugs, miconazole and clobetasol, increased the number of new oligodendrocytes, enhanced remyelination, and significantly reduced disease severity [82]. Treatment with human ESC‐derived oligodendrocytes functionally remyelinated the irradiated rat brain, as shown by complete recovery of cognitive deficits [83]. In transgenic AD mice, exercise prevented demyelination in the white matter, and improved their cognitive function [84]. Meanwhile, another study found oligodendrocyte induced by iPSCs transplantation ameliorated cognitive dysfunction in AD mice [85]. Most recently, enhancing myelin renewal can alleviate AD‐related cognitive impairment [86].

Targeted manipulations of cell‐type‐specific or tissue‐specific GRNs have allowed us to regenerate or restore cells lost to human injuries and degenerative diseases, such as regeneration of fat cells from myofibroblasts during wound healing and restoration of retinal neurons from Müller glia [87, 88]. GRNs control development via cell‐type‐specific gene expression and interactions between TFs and cis‐elements in the regulatory promoter regions. By integrating gene expression and chromatin accessibility data, we reconstructed gene regulatory networks for DEGs from PD cingulate cortex. Specifically, we built a sub‐network for myelination and gliogenesis that lay the foundation for remyelination by possible reactivation of its unique GRN.

A variety of animal models have been used to capture the clinical symptoms and investigate the aetiology of two most common neurodegenerative diseases: AD and PD. Current PD models can be divided into two categories: genetic and neurotoxic. None of these models recapitulates the age of onset, temporal progression, and the spectrum of problems and pathologies seen in the PD patients. It is expected that PD models would have a progressive loss of DA neurons. So far, many animal models have been generated and described for PD genes. The α‐synuclein transgenic mice, which expresses human mutant A53T α‐synuclein under the control of the mouse prion protein promoter, successfully developed α‐synuclein pathology and resembled many key features of progressive age‐dependent neurodegeneration as observed in PD patients [89, 90]. However, α‐synuclein based on transgenic mice exhibited little or no nigral degeneration. Similarly, for autosomal recessive genes, single knock‐out mice lacking PINK1, Parkin or DJ‐1 and even triple knock‐out mice lacking all these three recessive PD genes failed to show degeneration of DA neurons in the substantia nigra [91, 92, 93, 94]. Clearly, PD is a complex multifactorial disease, and a single model modality is not sufficient. During the past few decades, it is undeniable that studies and findings in PD animal models have played central roles for the development of new concepts and therapies used today. Even the challenges in PD modelling remain, we may have to accept the animal models are indispensable for future mechanistic studies and drug discoveries. Quite a few DEGs were validated in the two mouse models we used. Interestingly, some of them have been previously implicated in PD or neurodegeneration, which will be good candidates for further investigation of PD pathogenesis and aetiology using these PD models. It appeared that most DEGs involved in myelination and oligodendrocyte differentiation failed to be confirmed. We used 8‐month‐old Tg M83 mice for qPCRs. According to previous reports, Tg M83 mice developed motor phenotypes at 8 months of age [95]. However, non‐motor symptoms such as cognitive impairment and demyelination in the cingulate cortex for PD mouse model may occur in much later stages. Therefore, we may need to carefully select an appropriate animal model at elder ages when examining demyelination and myelin damage.

In conclusion, this study demonstrated the altered transcriptomic regulation of the cingulate cortex under PD condition that were not identified in previous studies. Our findings highlight how neurons and glia in the cingulate cortex respond to PD differently, as glutamate receptor signalling pathway upregulated and myelination network downregulated. Furthermore, integrative multi‐omics approaches were applied to construct GRNs that provide a transcriptional basis for remyelination therapies in the future.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Permission to reproduce materials from other sources

None.

Supplementary Material S1

Click here for additional data file.

Supplementary Material S2

Click here for additional data file.

Figure S1

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.