Transcriptomic signatures of brain regional vulnerability to Parkinson’s disease

The molecular mechanisms underlying the caudal-to-rostral progression of Lewy body pathology in Parkinson’s disease (PD) remain poorly understood. Here, we aimed to unravel transcriptomic signatures across brain regions involved in Braak Lewy body stages in non-neurological controls and PD donors. Using human postmortem brain datasets of non-neurological adults from the Allen Human Brain Atlas, we identified expression patterns related to PD progression, including genes found in PD genome-wide associations studies: SNCA, ZNF184, BAP1, SH3GL2, ELOVL7, and SCARB2. We confirmed these patterns in two datasets of non-neurological subjects (Genotype-Tissue Expression project and UK Brain Expression Consortium) and found altered patterns in two datasets of PD patients. Additionally, co-expression analysis across vulnerable regions identified two modules associated with dopamine synthesis, the motor and immune system, blood-oxygen transport, and contained microglial and endothelial cell markers, respectively. Alterations in genes underlying these region-specific functions may contribute to the selective regional vulnerability in PD brains.


6
BRGs present in UKBEC, 139 out of 317 (43.8%) negatively correlated BRGs and 400 out of 571 (70.1%) 119 positively correlated BRGs were differentially expressed between R1 and R6 (|FC R1-R6 | > 1, BH-corrected 120 P < 0.05, and Figure 2d). Second, we used RNA-sequencing (RNA-seq) data from 88-129 individuals in 121 the GTEx consortium 19  We next hypothesized that if the identified BRGs are associated with vulnerability to PD, they are also 128 indicative of vulnerability differences between PD patients and age-matched controls. To test this 129 R4-R6 also have higher expression in controls compared to PD patients. These findings thus support the 147 relation of BRGs with PD vulnerability encountered in brain regions of non-neurological individuals and 148 show how their expression may influence the vulnerability at a region-specific level as well as between 149 patients and controls. 150 as locomotory behavior and dopamine biosynthetic process, as well as diseases including 174 dopa−responsive dystonia, dystonialimb, and tremor, highlighting their role in motor circuitry. M47 was 175 enriched for endothelial cell markers and genes involved in immune response, blood coagulation, 176 interferon-gamma-mediated signaling pathway, and oxygen transport. This module was also enriched for 177 genes involved in auto-inflammatory or auto-immunity disorders, e.g. hypersensitivity, infection, and 178 inflammation. These modules and their associated pathways were associated with the preclinical stages 179 of PD, because of their higher expression in regions R1- R3. 180 Modules that were positively correlated with Braak stages were specifically enriched for neuronal markers 181 and related functions (e.g. axon, cell junction, and chemical synaptic transmission) reflecting higher 182 expression of these modules in the synapse-dense cerebral cortex. BRGs to be differentially expressed (BH-corrected P < 0.05 and β ≠ 0) between regions R1 and R6 after 194 correcting for five major cell-types (neurons, astrocytes, oligodendrocytes, microglia, and endothelial 195 cells). For example, the neuronal marker ADCY1 which was identified as a BRG remains differentially 196 expressed between regions R1 and R6 when corrected for neurons or other cell-types (Figure 4). 197 Similarly as for BRGs, PSEA analysis on all 23 Braak stage related co-expression modules showed 198 significant differential expression between regions R1 and R6 which cannot be fully explained by 199 differences in cellular composition. 200 In the PD datasets, not all BRGs were found significant after correction for cellular composition, however 201 smaller changes can be expected when comparing regions that are less distant (R1-R3 and R3-R4/R5). 202 Similar to the differential expression analysis without correction for cellular composition (Supplementary 203 Figure 3), PSEA revealed most changes between brain regions than between patients and controls 204 (Supplementary Figure 6). 205 Expression of several Parkinson's disease-implicated genes are related to 206 Braak staging 207 We found that the expression patterns of several PD-implicated genes, identified in the two most recent 208 PD genome-wide association studies 7,8 , were correlated with the Braak LB staging scheme. These 209 included BRGs (SCARB2, ELOVL7, SH3GL2, SNCA, BAP1, and ZNF184; Table 1  We further explored the relationship between SNCA expression and PD vulnerability in more detail. 213 observation suggests that lower SNCA expression indicates high vulnerability of brain regions to develop 217 LB pathology. We further validated this concept in two cohorts of PD patients in which SNCA expression 218 similarly increased across the medulla oblongata (R1), locus ceruleus (R2), and substantia nigra (R3) of 219 PD-, iLBD patients, and age-matched controls. SNCA was significantly lower expressed in region R1 220 compared to R2 and R3 in PD-and iLBD patients, but not in controls (Figure 5f). In the RNA-seq dataset, 221 SNCA was significantly lower expressed in the substantia nigra (R3) compared to the medial temporal 222 gyrus (R4/R5) in PD patients, but again not in controls ( Figure 5g). Altogether, SNCA expression patterns 223 could be replicated in brain regions of age-matched controls, however changes were larger between brain 224 regions in PD-and iLBD cases. We further assessed SNCA expression using PSEA in the AHBA ( Figure  225 5h) and found that changes were independent of neuronal or other cell-type densities when comparing 226 different brain regions. In PD datasets results were scattered and did not align between both PD 227 microarray and RNAseq datasets because of the small sample sizes and the comparison of different 228 brain regions (Supplementary Figure 8).

238
In PD, the progressive accumulation of LB pathology across the brain follows a characteristic pattern, 239 which starts in the brainstem and subsequently evolves to more rostral sites of the brain (Braak 240 ascending scheme) 1 . Using transcriptomic data of non-neurological brains, we identified genes (e.g. 241 SNCA, SCARB2, and ZNF184) and modules of co-expressed genes for which the expression decreased 242 or increased across brain regions defined by the Braak ascending scheme. Interestingly, these patterns 243 were disrupted in brains of patients with PD across regions that are preclinically involved in the 244 pathophysiology of PD. One gene co-expression module that showed higher expression in preclinically 245 involved regions was related to dopamine synthesis, locomotory behavior, and microglial and neuronal 246 activity. Another module was related to blood oxygen transport, the immune system, and may involve 247 endothelial cells. Our results highlight the complex genetic architecture of PD in which the combined 248 effects of genetic variants and co-expressed genes may underlie the selective regional vulnerability of the 249 brain. 250 Multiple studies suggests that a cytotoxic role and prion-like transfer of α-synuclein may contribute to its 251 progressive spread across the brain in PD, assuming a gain-of-function 3,27,28 . In line with this assumption 252 are reports of familial PD caused by SNCA multiplications, suggesting a SCNA dosage effect in causing 253 PD 5,29 . Interestingly, in contrast to the temporal and spatial pattern of the α-synuclein distribution 254 associated with the ascending Braak scheme in PD, the SNCA expression signature across brain regions 255 R1-R6 in non-neurological brains followed a reverse pattern with lowest expression in preclinically 256 involved regions (brainstem) and highest expression in clinically involved regions ( suppressor. Cancer-associated mutations within this gene were found to destabilize protein structure 282 promoting β-amyloid aggregation in vitro, which is the pathological hallmark in Alzheimer's disease 33 . 283 A number of functional pathways have been suggested to play a role in the pathogenesis of PD, such as 284 lysosomal function, immune system response, and neuroinflammation [6][7][8]34 . We identified modules of 285 genes that co-express across the six Braak stage-related regions and found they were enriched for genes 286 related to molecular processes that have been linked to the (pre)clinical symptoms and functional deficits 287 in PD. 288 One negatively correlated module M127 was enriched for genes related to functions and diseases 289 involving dopamine synthesis and motor functions. This module also contains the PD variant-associated 290 gene GCH1 (GTP cyclohydrolase 1) that is known to co-express with TH (tyrosine hydroxylase, the 291 enzyme responsible for converting tyrosine to L-3,4-dihydroxyphenylalanine (L-DOPA) in the dopamine 292 synthesis pathway) to enhance dopamine production and enable recovery of motor function in rat models 293 of PD 35 . In this study, both GCH1 and TH occur in M127 and thus were co-expressed across brain 294 regions involved in Braak stages supporting their interaction (Figure 6 and Supplementary Figure 9). The 295 higher expression in the more vulnerable brain regions R1-R3 indicates that GCH1, TH, and possibly 296 other genes within module M127 are essential to maintain dopamine synthesis that is affected in the early 297 Braak stages of PD. Indeed, by inhibiting TH activity, α-synuclein can act as a negative regulator of 298 dopamine release 26,36 . In this module, SLC18A2 (VMAT2; vesicular monoamine transporter 2) and 299 SLC6A3 (DAT; dopamine transporter) were also present, which are important for storage of dopamine 300 and transport in the cell 26 . Interestingly, dopamine may increase neuronal vulnerability, as was suggested 301 by an earlier study showing that α-synuclein is selectively toxic in dopaminergic neurons, and 302 neuroprotective in non-dopaminergic cortical neurons 37 . Cell-type marker enrichment showed that module 303 M127 was enriched for microglia-and neuronal markers, suggesting a role in neuroinflammation. α-304 Synuclein aggregates evoke microglia activation which in turn promotes aggregated protein propagation 305 to other brain regions, possibly even from the gut or periphery to the brain 27,34 . The higher expression of 306 microglial genes within module M127 may contribute to the higher vulnerability of brain regions affected 307 during preclinical stages to form protein aggregates. Further investigation of genes within module M127 will provide a better understanding of the molecular mechanisms underlying microglia activation, 309 dopaminergic pathways and motor functions. 310 Another negatively correlated module M47 was enriched for endothelial cell markers and genes involved 311 in functions and disorders that relate to the immune response and oxygen transport in blood. One 312 previous case-control study showed that anemia or low hemoglobin levels may precede the onset of 313 PD 38 . Several studies using blood transcriptomic meta-analysis revealed genes associated with 314 hemoglobin and iron metabolism were downregulated in PD patients compared to controls [39][40][41] . In our 315 study, several hemoglobin genes (HBD, HBB, HBA1, HBA2, and OASL) were also present in module 316 M47 of which HBD and HBB have been described to be highly interconnected with SNCA 41 . We also 317 found an association between the interferon-gamma-mediated signaling pathway and M47 in which OASL 318 also plays a role. Module M47 was negatively co-expressed with SNCA. Notably, a significant loss of 319 negative co-expression between SNCA and interferon-gamma genes in the substantia nigra has been 320 demonstrated in PD patients as compared to controls 42 . This loss may result from a downregulation of 321 genes within M47 in the substantia nigra of PD patients, similarly as was observed in PD blood 322 transcriptomics [39][40][41] . Therefore, these genes have the potential to serve as blood biomarkers for PD 323 vulnerability. Overall, these studies suggest that dysregulation of genes within module M47 involved in 324 blood-oxygen transport and the immune system influence brain regions to be selectively vulnerable to 325

PD. 326
Identification of transcriptomic features in regions or disease conditions may be confounded by changes 327 in cell-type composition. We used PSEA 24 to examine the impact of this confounding factor and found that 328 all 960 BRGs remained differentially expressed between regions R1 and R6 in AHBA. We also applied 329 PSEA in the two PD datasets that allowed us to examine cell-type specificity between regions as well as 330 between disease conditions. Although it is known that gene expression varies more between regions than 331 between disease conditions 22 , it is less clear how cell-type composition contributes to this variation. Here, 332 we found that regional comparisons yielded more significant results than when comparing disease changes were less dependent on cell-type abundance between regions than between patients and 335 controls. 336 In conclusion, we identified genes and pathways that may be important to maintain biological processes 337 in specific brain regions, but may also contribute to a higher selective vulnerability to PD. Our results 338 suggest that interactions between microglial genes and genes involved in dopamine synthesis and motor 339 functions, as well as between genes involved in blood-oxygen transport and the immune system may 340 contribute to the early involvement of specific brain regions in PD progression. Our observations highlight 341 a potential complex interplay of pathways in healthy brains and provide clues for future genetic targets 342 concerning the pathobiology in PD brains. 343

344
Methods and any associated references are available in the online version of the paper. 345 trans-synaptic alpha-synuclein transfer predicts regional atrophy in Parkinson disease.

486
The resulting modules of genes were subsequently analyzed to detect common biologically meaningful pathways.   Supplementary Table 5 and 7).

535
In the PD datasets, SNCA expression was tested for differential expression between (f, g) regions and conditions 536 (Supplementary Figure 6) and was considered significant when BH-corrected P < 0.05 (red).  Online methods 557 Allen Human Brain Atlas. To examine gene expression patterns across brain regions involved in 558 Parkinson's disease (PD), we used normalized gene expression data from the Allen Human Brain Atlas 559 (AHBA), a microarray data set of 3,702 anatomical brain regions from six non-neurological individuals (5 560 males and 1 female, mean age 42, range 24-57 years 15 ). We downloaded the data from 561 http://human.brain-map.org/. To filter and map probes to genes, the data was concatenated across the six 562 donors. We removed 10,521 probes with missing Entrez IDs, and 6,068 probes with low presence as they 563 were expressed above background in <1% of samples (PA-call containing presence/absence flag 15 ). The 564 remaining 44,072 probes were mapped to 20,017 genes with unique Entrez IDs using the collapseRows-565 function in R-package WGCNA v1.64.1 43 as follows: i) if there is one probe, that one probe is chosen, ii) if 566 there are two probes, the one with maximum variance across all samples is chosen 567 (method="maxRowVariance"), iii) if there are more than two probes, the probe with the highest 568 connectivity (summed adjacency) is chosen (connectivityBasedCollapsing=TRUE). Based on the 569 anatomical labels given in AHBA, samples were mapped to Braak stage-related regions R1-R6 as 570 defined by the BrainNet Europe protocol 18 and each region corresponds to one or multiple anatomical 571 structures. The locus ceruleus and pontine raphe nucleus are both part of the pontine tegmentum in R2. 572 Braak stage-related genes (BRGs). Two analysis methods were used to find genes for which the spatial 573 expression in AHBA is related to the spread of the disease: i) correlations between gene expression and 574 labels 1-6 according to their assignment to one of the Braak stage-related regions R1-R6, and ii) 575 differential expression between Braak stage-related regions R1 and R6. As the expression values were 576 log 2 -transformed, the mean difference between two regions was interpreted as the fold-change (FC). 577 Genes were assigned as BRGs based on the overlap of the top 10% (2,001) ranked genes with: i) highest 578 absolute correlation of gene expression and Braak stage labels, ii) highest absolute FC between R1 and To avoid capturing donor-specific changes, we applied the correlation and differential expression 581 analyses for each of the six brain donors separately, and effect sizes were then combined by meta-582 analysis (metafor R-package 2.0). A random effects model was applied which assumes that each brain is 583 considered to be from a larger population of brains and therefore takes the within-brain and between-584 brain variance into account. The between-brain variance (tau 2 ) was estimated with the Dersimonian-585 Delaird model. Variances and confidence intervals were obtained using the escalc-function. Correlations 586 were Fisher-transformed (z = arctanh(r)) to obtain summary estimates, which were then back-transformed 587 to correlation values ranging between -1 and +1. P-values were BH-corrected for all 20,017 genes. The 588 significance of summary effect sizes (correlations and FCs) was assessed through a two-sided t-test (H 0 : 589 FC=0; unequal variances). The weights used in the meta-analysis are based on the non-pooled 590 expression variance in R1-R6. 591 UK Brain Expression Consortium (UKBEC). UKBEC 20 (http://www.braineac.org) contains microarray 592 expression data from 10 brain regions of 134 non-neurological donors (74.5% males, mean age 59, range 593 16-102 years) for which their control status was confirmed by histology. We used the biomaRt R-package 594 version 2.38 44 to map Affymetrix probe IDs from UKBEC to gene Entrez IDs; 262,134 out of 318,197 595 probes could be mapped. Similar as for AHBA, expression data for all probes and samples was 596 concatenated across the 10 brain regions before mapping probes to 18,333 genes with unique Entrez IDs 597 using the collapseRows-function.