https://orcid.org/0000-0002-8818-5382
Abstract
Extensive studies have demonstrated significant gray matter atrophy in patients with Parkinson’s disease (PD); however, the underlying gene expression mechanisms remain largely unknown. To comprehensively characterize the gray matter volume alterations in PD patients, we conducted a neuroimaging meta-analysis and validated the observed atrophic phenotypes in an independent dataset. Leveraging the Allen Human Brain Atlas (AHBA), we linked brain transcriptomic data to neuroimaging phenotypes to identify genes associated with PD-related gray matter atrophy. Further enrichment analyses and functional characterization explored the potential roles of these correlated genes in disease pathology. Both the neuroimaging meta-analysis and independent dataset analysis consistently revealed significant gray matter atrophy in PD, particularly in the superior temporal gyrus, highly associated with sensory and motor functions. Spatial transcriptome-neuroimaging correlation analysis identified 1,952 overlapping genes whose expression levels were significantly correlated with the spatial distribution of gray matter atrophy in PD patients. These genes were enriched in several key biological processes and molecular pathways, exhibiting region- and cell type-specific expression, particularly in dopaminergic receptor neurons of brain tissue. This study delineates the spatial distribution of gray matter atrophy in PD and suggests that this neurodegenerative phenotype may result from complex interactions among multiple functionally relevant genes.
Introduction
Parkinson’s disease (PD) is a fast-growing neurodegenerative condition on a global scale, thus giving rise to a considerable societal burden (Bloem et al. 2021). Notably, PD represents clinically, pathologically, and genetically heterogeneous, which resists distillation into a single, cohesive disorder (Ye et al. 2023). Similarly, its pathology involves extensive regions of the nervous system, various neurotransmitters, and protein aggregates other than just Lewy bodies (Kalia and Lang 2015). As PD could be a fragmented collection of many different conditions with variable clinical and pathologic overlap, it may be more beneficial to gain insights into the pathogenesis from both clinical and genetic perspectives, as well as molecular pathology.
Evidence from both neuroimaging and clinicopathologic studies strongly suggests that the loss of synapses and axonal degeneration contribute to initial PD motor manifestations, even preceding neuronal loss (O'Keeffe and Sullivan 2018). Much effort has been devoted to exploring the exact morphological changes in the brain, resulting in heterogeneous findings and controversial concepts (Jellinger 2012). Thereinto, voxel-based morphometry (VBM) is a powerful tool for identifying brain structural abnormalities in vivo. With respect to PD, numerous relevant morphological studies using VBM methodology have been conducted to examine gray matter volume (GMV) abnormalities between PD patients and healthy people (Burton et al. 2004; Nagano-Saito et al. 2005; Sanchez-Castaneda et al. 2009; Kostić et al. 2010; Meppelink et al. 2011; Archibald et al. 2013; Lin et al. 2015; Chen et al. 2019). However, the present results are mixed and conflicting across studies, along with the lack of in-depth interpretation and exploration.
Genetics is an undeniable factor when discussing the causes of PD (Bloem et al. 2021). We have now recognized that virtually all PD cases likely have detectable genetic influences, with the specific genetic variants involved in individual cases varying by frequency and effect size. Beginning with the landmark 1997 discovery of mutations in synuclein-alpha (SNCA) (Polymeropoulos et al. 1997), the last 27 years have witnessed remarkable progress, with approximately 100 distinct genes or loci having now been definitively linked to PD susceptibility (Blauwendraat et al. 2020). Moreover, recent genome-wide genotyping in large cohorts has also identified 90 independent genome-wide variants and enabled the estimation of PD heritability, with up to 36% of overall disease risk attributed to common genetic variation (Nalls et al. 2019). Altogether, these studies have substantially expanded the resources available for future investigations into potential PD pathology and intervention.
Taken together, it is well established that brain structural damage is a typical and stable neuropathological feature of PD. Nevertheless, the genetic mechanisms underlying this neurobiological phenotype are far from being understood. A combined analysis of brain imaging data and brain-wide gene expression atlases such as the Allen Human Brain Atlas (AHBA) has given rise to the emergent domain of neuroimaging transcriptomics (Zarkali et al. 2020; Keo et al. 2021; Thomas et al. 2021). This approach provides a feasible route toward the identification of genes with spatial profiles of regional expression that track anatomical variations in a certain neuroimaging phenotype. Furthermore, the stereotaxic mapping of these brain tissue samples has made it feasible to test the associations between spatial patterns in gene expression and spatial profiles in neuroimaging phenotypes given the common coordinate system (Arnatkeviciute et al. 2019), contributing to the rapidly growing field of neuroimaging transcriptomics that has attempted to bridge the gap between microscale molecular function and macroscale brain architecture (Fang et al. 2023; Sun et al. 2023).
The objective of our study was to investigate the genetic mechanisms underlying GMV alterations in PD. To achieve this, we first characterized GMV changes through a neuroimaging meta-analysis and a morphological study using the VBM method in an independent dataset, ensuring a balance between generalizability and data authenticity. Furthermore, we combined the AHBA to perform a transcriptome-neuroimaging spatial association analysis to identify genes whose expression levels were associated with gray matter atrophy in PD. Finally, an array of post hoc analyses (e.g., functional enrichment, specific expression, protein-protein interaction (PPI), and behavioral relevance analyses) were conducted to investigate the functional features of these identified genes. A schematic overview of the study design and analysis pipeline is shown in Fig. 1.
A schematic overview of the study design and analysis pipeline. Abbreviations: GM, gray matter; HC, healthy control; PD, Parkinson’s disease.
Methods
Literature search and data extraction
Our meta-analysis was carried out strictly following the Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) guidelines to ensure comprehensive and transparent reporting of methods and results (Appendix in Supplementary Materials) (Moher et al. 2009). Two researchers (Y.J. and M.X.) independently conducted a systematic search to determine the related studies in PubMed and Web of Science before June 26, 2024. We used the following keywords to search for suitable studies: (“voxel-based morphometry” OR “VBM” OR “morphometry” OR “gray matter” OR “grey matter”) AND (“Parkinson’s disease” OR (“Parkinson” AND “disease”) OR “Parkinson disease” OR “PD”). Moreover, the meta-analysis was preregistered on PROSPERO (https://www.crd.york.ac.uk/PROSPERO/, NO: CRD42023439786).
Study screening and data extraction were performed separately by two authors (Y.J. and M.X.). Any disagreement was resolved by the adjudicating senior investigator (F.X.). All studies were included according to the following criteria: (i) studies explored GMV changes in patients with PD; (ii) patients diagnosed as PD according to the criteria; (iii) comparisons were performed between PD patients and controls; (iv) voxel-wise comparisons were carried out within whole-brain gray matter rather than regions of interest (ROIs); (v) statistical results were reported in a stereotactic space (Montreal Neurological Institute (MNI) or Talairach coordinates). The exclusion criteria were as follows: (i) studies were not original papers (e.g., editorials, letters to the editor, review articles, and meta-analyses); (ii) studies were case reports; (iii) studies involved animal experiments; and (iv) data could not be obtained from the published investigations or after contacting the authors. We also excluded studies whose data overlapped those of other studies (e.g., from the same institution or authors), with the most recent or complete report included in our meta-analysis. In addition, we manually searched the reference lists of the selected articles and PD-related meta-analyses for additional qualifying studies.
For each included study, we recorded the following information: first author, sample size, sex, mean age, education, illness duration, UPDRS-III (Unified Parkinson’s Disease Rating Scale-Part III) scores, software used for VBM, and magnetic resonance imaging scanner. The research was conducted strictly in line with PRISMA guidelines (Moher et al. 2009).
Neuroimaging meta-analysis of GMV differences
A voxel-wise meta-analysis of GMV differences between PD patients and healthy controls (HCs) was conducted using Seed-based d Mapping with Permutation of Subject Images (SDM-PSI, version 6.22). In contrast to the conventional activation likelihood estimation approach, SDM-PSI combines standard effect size and variance-based meta-analytic computations to permit researchers to combine peak coordinates and statistical parametric maps (Albajes et al. 2019). For each study, we first extracted the peak coordinates and corresponding effect sizes (such as t values) of clusters with reported GMV differences between PD patients and HC, z or P-values were then converted to t values using SDM online conversion utilities (https://www.sdmproject.com/utilities/?show=Statistics). Then, using a Gaussian kernel of 20 mm full width at half maximum, a standard MNI map of GMV differences between PD patients and HCs was created independently for each research (Radua et al. 2012). The final map of GMV differences (z map) between groups for all included studies was created by combining these maps using a standard random-effects model accounting for sample size, intrastudy variability, and between-study heterogeneity. Multiple comparisons were corrected using the cluster-level family-wise error (FWE) approach, yielding a corrected cluster-level significance of P < 0.05 (Lieberman and Cunningham 2009; Radua et al. 2012).
Several supplementary analyses were pursued to test the robustness and reliability of our meta-analysis results. Potential publication bias was assessed by two tests for small-study effect and excess significance, respectively. Finally, the I2 statistic was calculated to describe which percentage of the variability between studies might be due to between-study heterogeneity, with I2 = 25%, 50%, and 75% indicating low, medium, and high heterogeneity, respectively (Higgins and Thompson 2002). Furthermore, we have conducted meta-regression analyses to examine the confounders, including the UPDRS-III total scores and the duration of the disease, which may have affected the results.
Participants recruited in the independent dataset
Group GMV differences were also investigated in an independent dataset of 48 PD patients and 26 HCs. Prospective and continuous recruitment of PD patients attending the PD Clinic at the Affiliated Wuxi People’s Hospital of Nanjing Medical University was conducted between August 2021 and April 2023. The diagnosis of PD was based on the clinical diagnostic criteria for PD (Postuma et al. 2015) by senior neurologists. A senior neurologist (L.Z.) initially diagnosed all subjects, and their evaluations were further assessed by a neurologist professor (F.W.). PD patients were included if they (i) were coincident with the PD diagnosis standard according to the clinical diagnostic criteria for PD (Postuma et al. 2015); (ii) were aged 40-80 years old; (iii) were right-handers; (iv) were able to perform MRI scans and finish the neurological and psychological assessment; and (v) were willing to take part in this research. PD patients were extra excluded: (i) inability to cooperate with clinical assessment or MRI examination; (ii) maximum head motion ≥2.5 mm or 2.5°; and (iii) history of severe neurological or cerebrovascular diseases besides PD. All dopaminergic therapy was withheld for at least 12 h before the MRI scan to alleviate the impact of drugs. Healthy subjects matched for age and sex were enrolled as controls. All participants gave their written informed consent in accordance with the Declaration of Helsinki after receiving approval (KY21133) from the ethics committee of the Affiliated Wuxi People’s Hospital of Nanjing Medical University.
Magnetic resonance images were acquired using a 3.0T MRI scanner (Magnetom 3T Siemens, Prisma, Germany). We obtained T1-weighted anatomical data using a volumetric 3D-magnetization prepared rapid acquisition gradient echo (3DT1WI MP-RAGE) sequence, with rigorous quality control throughout the process. When shifting lying positions, we used a foam pad to minimize head motion and earplugs to reduce scanner noise. In addition, all the participants were instructed to stay awake, close their eyes, and try not to think about anything during the examination (Fox et al. 2005). The detailed scanning parameters were as follows: TR = 2,300 ms, TE = 2.98 ms, TI = 900 ms, FA = 9°, slice thickness = 1 mm, slices per slab = 192, FOV = 256 × 256 mm2, matrix size = 256 × 256, and voxel size = 1 × 1 × 1 mm3. The 3D-T1WI MP-RAGE scanning duration was 5 min and 30 s.
GMV analysis in the independent dataset
We applied the Statistical Parametric Mapping program (SPM12) to run the VBM analysis on the MATLAB R2018b platform to explore structural differences among two groups in the independent dataset (Kurth et al. 2015). The VBM analysis consists of the following three steps: (i) segmenting the 3D-T1WI MP-RAGE data of each participant into the cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM); (ii) spatially normalizing and transforming the GM images into MNI space; and (iii) spatial smoothing at 8 × 8 × 8 mm3 FWHM Gaussian kernel. Total intracranial volumes (TIVs) were calculated by adding up the total voxels of GM, WM, and CSF in the native space separately. For the independent dataset, GMV differences between PD patients and HC were investigated using a two-sample t-test in a voxel-wise manner, with age, sex, and TIVs as covariates, resulting in a statistical t map. Multiple comparisons were corrected using the FWE approach, with a cluster-defining threshold of P < 0.001 and a corrected cluster-level significance of P < 0.05.
Brain gene expression data processing
Brain gene expression data were extracted from the open AHBA dataset (http://www.brain-map.org) (Hawrylycz et al. 2012; Hawrylycz et al. 2015). The AHBA has collected expression data of over 20,000 genes in 3,702 spatially distinct tissue samples from six donated human brains (Hawrylycz et al. 2012). The dataset was collected from six human post-mortem donors (Supplementary Table 1). The original expression data of over 20,000 genes via 3,702 spatially separate brain tissue samples was analyzed in accordance with a newly proposed pipeline (Arnatkeviciute et al. 2019). Using the Re-Annotator package, we particularly first updated the most probe-to-gene annotations based on the most up-to-date information from the National Center for Biotechnology Information (Arloth et al. 2015). We eliminated probes that did not outperform the background noise in at least 50% of samples from all donors using intensity-based filtering. We further applied RNA-seq data as a reference to pick probes since multiple probes were used to quantify the expression level of a single gene. We calculated the correlations between microarray and RNA-seq expression measures for the remaining genes after removing genes whose expression patterns do not overlap across RNA-seq and microarray datasets. After eliminating probes with low correlations (r < 0.2), a representative probe for a gene was chosen on the basis of the strongest correlation to the RNA-seq data. This study only used tissue samples from the left cerebral cortex. For one thing, all six donors had expression data in the left hemisphere, while only two had samples in the right. Another issue is that the inclusion of subcortical samples may generate biases due to the large disparities in gene expression between cortical and subcortical regions (Hawrylycz et al. 2012). Following the aforementioned processing procedure, we obtained normalized expression data for 5,013 genes across 1,280 tissue samples by applying the standard gray matter mask of the left cerebral cortex (Supplementary Fig. 1). Subsequently, we restricted all analyses—including the neuroimaging meta-analysis, VBM analysis, group-level comparison, and transcriptome-neuroimaging analysis—to the samples within this mask, resulting in a final sample × gene matrix of 773 × 5,013.
To account for both potential between-sample differences and donor-specific effects in gene expression, we applied the scaled robust sigmoid normalization approach to perform within-sample cross-gene and within-gene cross-sample normalization. Differential stability (DS) is an indicator that measures the consistency of regional variance among donor brains. Previous studies have found that genes with high DS scores exhibit more consistent spatial expression patterns among donors (Hawrylycz et al. 2015). We only picked genes with relatively more conserved expression patterns for analysis since gene expression conservation across subjects is necessary for the transcription–neuroimaging spatial correlations. In order to manage, we ranked the genes according to their DS values and selected the top 50% with the highest DS values for the primary analysis. Furthermore, we performed sensitivity analyses using two additional DS cutoff thresholds (top 40% and 60%) to assess the impact of various DS threshold selections.
Transcription-neuroimaging association analysis
To delineate the group GMV difference in a specific brain tissue sample, we created a spherical region (radius = 3 mm) centered on the sample’s MNI coordinate and extracted the average z value of voxels within the sphere from the meta-analysis z map. Then, gene-wise cross-sample Pearson’s correlations (773 samples) between gene expression and z values were carried out, resulting in 5,013 correlation coefficients. Multiple comparisons were corrected using the Benjamini and Hochberg method for false discovery rate (FDR-BH) (P < 0.05). The t map in the independent dataset was also subjected to the same analysis mentioned above. Then, we restricted our analysis to the genes whose expression levels were related to changes in GMV in PD patients and had significant spatial correlations with both the meta-analysis z map and the independent dataset t map.
Spatially constrained permutation test
In order to test the statistical significance of our findings, a spatially constrained permutation was conducted to determine whether the number of our discovered genes was significantly higher than the random level. The standard nonparametric null (e.g., randomly shuffling the sample labels) is strongly violated by the spatial autocorrelation of brain maps given that transcriptional data are spatially autocorrelated; for example, nearby anatomical regions tend to have more similar patterns of gene expression than spatially distant regions, resulting in greater family-wise error rates (Markello and Misic 2021). To deal with this issue, we conducted the permutation test using a spatially constrained null model suggested by Burt et al. (2020). This method is implemented in an open-access, Python-based software package, BrainSMASH: Brain Surrogate Maps with Autocorrelated Spatial Heterogeneity (https://github.com/murraylab/brainsmash). It can simulate volumetric surrogate brain maps, which utilize 3D Euclidean distance between regions to preserve spatial autocorrelation. We applied this approach to produce spatial autocorrelation-preserving surrogate maps for each gene to correct the spatial autocorrelation in transcriptional data. Using the same procedure as described above, these surrogate maps were utilized to re-identify genes linked to GMV variations in PD patients. Then, the number of genes discovered in each test was recorded as we performed this process 5,000 times to create a null distribution. Finally, to test whether our results varied from randomness, we compared the number of genes obtained through the real data with this null distribution.
Gene enrichment, functional feature and protein-protein interaction analysis
A series of enrichment analyses were performed for the identified genes associated with GMV differences in PD patients. First, functional annotation was performed using the ToppGene portal (https://toppgene.cchmc.org/) (Chen et al. 2009). The biological functions, including molecular functions (MFs), biological processes (BPs), cellular components (CCs), and pathways (PWs) were identified using the gene ontology (GO) and pathway databases. Second, to identify the specific tissues, cell types, and developmental stages in which these genes were overrepresented, we used online tissue-specific (http://genetics.wustl.edu/jdlab/tsea/) and cell type-specific expression analysis tools (http://genetics.wustl.edu/jdlab/csea-tool-2/) (Dougherty et al. 2010). To estimate the chance that a gene was specifically expressed, a specificity index probability (pSI) with four pSI thresholds (0.05, 0.01, 0.001, and 0.0001) was utilized (Xu et al. 2014). Finally, using 20,585 genes with distinct Entrez IDs in the AHBA as the background list, we examined the overlap between the genes associated with GMV changes in PD patients discovered in the current study and 90 independent genome-wide significant association signals identified through the largest genetic meta-analysis (Nalls et al. 2019). Fisher’s exact tests were employed to determine the statistical significance of the aforementioned enrichment analyses. With a corrected P-value (q) of 0.05, multiple testing was corrected using FDR-BH correction.
PPI analysis was carried out using STRING v11.0 to find out whether these identified genes could build a PPI network with the highest confidence interaction score of 0.9 (Szklarczyk et al. 2023). Hub genes were classified as genes with the top 10% of the most significant degree values (e.g., the number of edges connecting two genes). Moreover, the hub genes with the highest degree values had their spatial-temporal expression trajectory defined using the Human Brain Transcriptome database (http://hbatlas.org/).
Behavioral relevance analysis
We investigated the associations of these identified genes related to GMV changes in PD with behavioral domains from Neurosynth (https://neurosynth.org/), a well-validated and publicly available platform for meta-analysis of neuroimaging literature, to capture the behavioral relevance of the identified genes connected to GMV variations in PD patients (Yarkoni et al. 2011). Nearly all aspects of human behavior are covered by the activation maps of 1,335 behavioral terms contained in the Neurosynth database. Cross-sample correlation analysis between activation levels and gene expression measures was performed for each term, yielding a set of correlation coefficients corresponding to the genes. We averaged the absolute values of these correlation coefficients (|r| mean) to measure the degree to which these genes were associated with each term, given that high positive and negative correlation coefficients both show that a gene contributes to a behavioral term. Finally, the behavioral terms were ranked by |r| mean, and the terms with the highest |r| mean were chosen to represent the behavioral relevance of the genes associated with GMV changes in PD patients. For visualization and interpretation, a criterion of |r| mean > 0.2 was employed in this case. Using the same analytic pipeline applied to the Neurosynth, we also tested the associations of gene expression with 63 behavioral domains via the BrainMap (http://brainmap.org/) (Lancaster et al. 2012).
Results
Gray matter atrophy in PD from meta-analysis
Following the extensive literature review and selection process, 1,831 PD patients and 1,378 HCs from 44 studies were included in our neuroimaging meta-analysis (Supplementary Fig. 2). The demographic and clinical characteristics of these participants are listed in Supplementary Table 2. In comparison to HCs, PD patients showed prominent gray matter atrophy in the bilateral superior temporal gyrus and inferior parietal gyrus, the right angular gyrus, and the left superior frontal gyrus (P < 0.05, FWE-corrected) (Table 1 and Fig. 2a). Furthermore, publication bias analysis (including the small-study effect test, excess significance test, and I2 statistic) demonstrated the reliability and validity of our meta-analysis results (Supplementary Table 3 and Supplementary Fig. 3).
Brain regions with gray matter atrophy between PD patients and HCs identified from the neuroimaging meta-analysis a) and the independent dataset b). The color bar represents the level of gray matter atrophy, respectively. The differential brain maps were both visualized with the BrainNet Viewer (Xia et al. 2013). Abbreviations: HC, healthy controls; L, left; PD, Parkinson’s disease; R, right.
Brain regions showing gray matter atrophy between PD patients and HCs identified in the neuroimaging meta-analysis.
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | SDM-z | P-value |
|---|---|---|---|---|
| R-superior temporal gyrus | 7815 | 50, −2, −4 | −9.525 | 1.00 × 10−3 |
| R-angular gyrus | 1427 | 46, −66, 38 | −6.842 | 1.00 × 10−3 |
| R-inferior parietal | 585 | 52, −36, 46 | −5.923 | 2.00 × 10−3 |
| L-superior temporal gyrus | 404 | −48, 12, −14 | −6.582 | 1.00 × 10−3 |
| L-inferior parietal | 236 | −48, −28, 46 | −5.28 | 1.50 × 10−2 |
| L-superior frontal gyrus | 210 | 2, 56, 16 | −5.771 | 1.90 × 10−2 |
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | SDM-z | P-value |
|---|---|---|---|---|
| R-superior temporal gyrus | 7815 | 50, −2, −4 | −9.525 | 1.00 × 10−3 |
| R-angular gyrus | 1427 | 46, −66, 38 | −6.842 | 1.00 × 10−3 |
| R-inferior parietal | 585 | 52, −36, 46 | −5.923 | 2.00 × 10−3 |
| L-superior temporal gyrus | 404 | −48, 12, −14 | −6.582 | 1.00 × 10−3 |
| L-inferior parietal | 236 | −48, −28, 46 | −5.28 | 1.50 × 10−2 |
| L-superior frontal gyrus | 210 | 2, 56, 16 | −5.771 | 1.90 × 10−2 |
Abbreviations: HC, healthy control; L, left; MNI, Montreal Neurological Institute; PD, Parkinson’s disease; R, right; SDM, Seed-based d-Mapping.
Brain regions showing gray matter atrophy between PD patients and HCs identified in the neuroimaging meta-analysis.
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | SDM-z | P-value |
|---|---|---|---|---|
| R-superior temporal gyrus | 7815 | 50, −2, −4 | −9.525 | 1.00 × 10−3 |
| R-angular gyrus | 1427 | 46, −66, 38 | −6.842 | 1.00 × 10−3 |
| R-inferior parietal | 585 | 52, −36, 46 | −5.923 | 2.00 × 10−3 |
| L-superior temporal gyrus | 404 | −48, 12, −14 | −6.582 | 1.00 × 10−3 |
| L-inferior parietal | 236 | −48, −28, 46 | −5.28 | 1.50 × 10−2 |
| L-superior frontal gyrus | 210 | 2, 56, 16 | −5.771 | 1.90 × 10−2 |
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | SDM-z | P-value |
|---|---|---|---|---|
| R-superior temporal gyrus | 7815 | 50, −2, −4 | −9.525 | 1.00 × 10−3 |
| R-angular gyrus | 1427 | 46, −66, 38 | −6.842 | 1.00 × 10−3 |
| R-inferior parietal | 585 | 52, −36, 46 | −5.923 | 2.00 × 10−3 |
| L-superior temporal gyrus | 404 | −48, 12, −14 | −6.582 | 1.00 × 10−3 |
| L-inferior parietal | 236 | −48, −28, 46 | −5.28 | 1.50 × 10−2 |
| L-superior frontal gyrus | 210 | 2, 56, 16 | −5.771 | 1.90 × 10−2 |
Abbreviations: HC, healthy control; L, left; MNI, Montreal Neurological Institute; PD, Parkinson’s disease; R, right; SDM, Seed-based d-Mapping.
After careful processing of the clinical indicators, the meta-regression analyses indicated that those patients with longer disease durations tended to have more apparent gray matter atrophy in the left superior medial frontal gyrus, left fusiform gyrus, bilateral middle cingulate gyrus, and left supplementary motor area. Higher UPDRS-III scores were associated with gray matter atrophy in the left insula, bilateral middle temporal gyrus, and right superior and inferior temporal gyri. However, these findings should be interpreted cautiously because no clusters survived in the correlation brain maps when a more rigorous correction (50 times multiple permutations) was applied.
Gray matter atrophy in PD from the independent dataset
The demographic and clinical characteristics of the participants in the independent dataset are presented in Table 2. There were no significant differences in age, sex, or year of education between PD patients and HC. In patients with PD, the voxel-wise two-sample t-test showed merely gray matter atrophy, comparable to the results in the meta-analysis (Fig. 2b and Table 3). Remarkably, PD patients exhibited gray matter atrophy in the left superior temporal gyrus and amygdala, partially consistent with findings from the meta-analysis (P < 0.05, cluster-level FWE-corrected).
Demographic and clinical characteristics of the independent dataset.
| Characteristics | PD | HC | Statistics | P-value |
|---|---|---|---|---|
| Number of subjects | 46 | 28 | ||
| Sex (M/F) | 24/22 | 15/13 | X2 = 0.014 | 0.907 |
| Age (years) | 64.18 ± 9.38 | 64.18 ± 8.14 | t = 0.002 | 0.998 |
| Education (years) | 10.23 ± 2.88 | 9.81 ± 3.39 | t = 0.569 | 0.571 |
| Disease duration (years) | 4.31 ± 3.05 | NA | ||
| H&Y | 2.09 ± 0.78 | NA | ||
| UPDRS-III scores | 23.77 ± 12.06 | NA |
| Characteristics | PD | HC | Statistics | P-value |
|---|---|---|---|---|
| Number of subjects | 46 | 28 | ||
| Sex (M/F) | 24/22 | 15/13 | X2 = 0.014 | 0.907 |
| Age (years) | 64.18 ± 9.38 | 64.18 ± 8.14 | t = 0.002 | 0.998 |
| Education (years) | 10.23 ± 2.88 | 9.81 ± 3.39 | t = 0.569 | 0.571 |
| Disease duration (years) | 4.31 ± 3.05 | NA | ||
| H&Y | 2.09 ± 0.78 | NA | ||
| UPDRS-III scores | 23.77 ± 12.06 | NA |
Abbreviations: HC, healthy controls; H&Y, Hoehn & Yahr scales; NA, not applicable; PD, patients with Parkinson’s disease; UPDRS-III, Unified Parkinson’s Disease Rating Scale.
Demographic and clinical characteristics of the independent dataset.
| Characteristics | PD | HC | Statistics | P-value |
|---|---|---|---|---|
| Number of subjects | 46 | 28 | ||
| Sex (M/F) | 24/22 | 15/13 | X2 = 0.014 | 0.907 |
| Age (years) | 64.18 ± 9.38 | 64.18 ± 8.14 | t = 0.002 | 0.998 |
| Education (years) | 10.23 ± 2.88 | 9.81 ± 3.39 | t = 0.569 | 0.571 |
| Disease duration (years) | 4.31 ± 3.05 | NA | ||
| H&Y | 2.09 ± 0.78 | NA | ||
| UPDRS-III scores | 23.77 ± 12.06 | NA |
| Characteristics | PD | HC | Statistics | P-value |
|---|---|---|---|---|
| Number of subjects | 46 | 28 | ||
| Sex (M/F) | 24/22 | 15/13 | X2 = 0.014 | 0.907 |
| Age (years) | 64.18 ± 9.38 | 64.18 ± 8.14 | t = 0.002 | 0.998 |
| Education (years) | 10.23 ± 2.88 | 9.81 ± 3.39 | t = 0.569 | 0.571 |
| Disease duration (years) | 4.31 ± 3.05 | NA | ||
| H&Y | 2.09 ± 0.78 | NA | ||
| UPDRS-III scores | 23.77 ± 12.06 | NA |
Abbreviations: HC, healthy controls; H&Y, Hoehn & Yahr scales; NA, not applicable; PD, patients with Parkinson’s disease; UPDRS-III, Unified Parkinson’s Disease Rating Scale.
Brain regions showing gray matter atrophy between PD patients and HCs identified from the independent dataset.
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | t value |
|---|---|---|---|
| L-superior temporal gyrus | 646 | −58.5, −19.5, 10.5 | 5.1296 |
| L-middle temporal gyrus | 119 | ||
| L-amygdala | 267 | −16.5, 0, −16.5 | 4.1336 |
| L-parahippocampal gyrus | 155 | ||
| L-hippocampus | 154 | ||
| L-fusiform | 128 | ||
| L-middle temporal gyrus | 100 | ||
| L-inferior temporal gyrus | 97 | ||
| L-superior temporal gyrus | 89 |
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | t value |
|---|---|---|---|
| L-superior temporal gyrus | 646 | −58.5, −19.5, 10.5 | 5.1296 |
| L-middle temporal gyrus | 119 | ||
| L-amygdala | 267 | −16.5, 0, −16.5 | 4.1336 |
| L-parahippocampal gyrus | 155 | ||
| L-hippocampus | 154 | ||
| L-fusiform | 128 | ||
| L-middle temporal gyrus | 100 | ||
| L-inferior temporal gyrus | 97 | ||
| L-superior temporal gyrus | 89 |
Abbreviations: HC, healthy control; L, left; MNI, Montreal Neurological Institute; PD, Parkinson’s disease.
Brain regions showing gray matter atrophy between PD patients and HCs identified from the independent dataset.
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | t value |
|---|---|---|---|
| L-superior temporal gyrus | 646 | −58.5, −19.5, 10.5 | 5.1296 |
| L-middle temporal gyrus | 119 | ||
| L-amygdala | 267 | −16.5, 0, −16.5 | 4.1336 |
| L-parahippocampal gyrus | 155 | ||
| L-hippocampus | 154 | ||
| L-fusiform | 128 | ||
| L-middle temporal gyrus | 100 | ||
| L-inferior temporal gyrus | 97 | ||
| L-superior temporal gyrus | 89 |
| Brain regions | Cluster size (voxels) | Peak MNI coordinate | t value |
|---|---|---|---|
| L-superior temporal gyrus | 646 | −58.5, −19.5, 10.5 | 5.1296 |
| L-middle temporal gyrus | 119 | ||
| L-amygdala | 267 | −16.5, 0, −16.5 | 4.1336 |
| L-parahippocampal gyrus | 155 | ||
| L-hippocampus | 154 | ||
| L-fusiform | 128 | ||
| L-middle temporal gyrus | 100 | ||
| L-inferior temporal gyrus | 97 | ||
| L-superior temporal gyrus | 89 |
Abbreviations: HC, healthy control; L, left; MNI, Montreal Neurological Institute; PD, Parkinson’s disease.
Genes associated with gray matter atrophy in PD
Based on transcriptome-neuroimaging spatial correlation analyses, we found that the meta-analysis z map and the independent dataset t map were associated with expression measures of 4,287 and 2,168 genes, respectively (P < 0.05, FDR-corrected), with 1,952 overlap genes (Supplementary Dataset 1). Further analyses were conducted using these identified overlap genes, of which expression levels were thought to be closely associated with gray matter atrophy in PD patients. The rigorous spatially constrained permutation test revealed that none of the 5,000 permutations produced more genes than those found using the real data (Pperm < 0.0002), indicating that our results were by no means random (Supplementary Fig. 4). Additionally, we observed similar results using two other DS thresholds of 40% (overlap ratio: 87.52%) and 60% (overlap ratio: 93.12%) (Supplementary Table 4 and Supplementary Dataset 2). The 90 independent genome-wide significant risk signals and the identified 1,952 genes associated with gray matter atrophy in PD patients were significantly overlapping by use of a Fisher’s exact test (14 overlap genes, odds ratio = 1.78, P = 0.0405). Specifically, these 14 overlapping genes were listed below: CAMK2D; CNTN1; DDRGK1; ELOVL7; FAM49B; GPNMB; LRRK2; MBNL2; RIT2; SATB1; SEMA4A; SNCA; SYT17; WNT3.
Gene functional enrichment
We conducted functional enrichment analysis utilizing the ToppGene portal to clarify the biological processes and pathways of the identified genes involved in gray matter atrophy in PD patients. The results of functional enrichment are illustrated in Supplementary Fig. 5 along with Supplementary Dataset 3. With regard to GO, these genes were enriched for MFs including gated and voltage-gated channel activity, cytoskeletal protein binding, channel regulator and inhibitor activity, calcium, sodium and potassium channel activity, and calcium ion binding; for BPs including the central nervous system, neuron, brain and sensory organ development, cell morphogenesis involved in neuron differentiation, regulation of nervous system development, positive regulation of signal transduction, developmental process, and neurogenesis and axonogenesis and, for CCs, including neuron projection, synapse, axon, dendrite, neuron to neuron synapse, neuronal cell body, GABA-ergic synapse, neuro projection cytoplasm, spine, and projection membrane. As to pathways, these genes were found to be apparently enriched for the neuronal system, potassium channels, cell-cell communication, transmission across chemical synapses, neurotoxicity of clostridium toxins, neurotransmitter receptors, and transmission across chemical synapses.
Tissue-, cell type-, and temporal-specific expressions
The specific expression results of the 1,952 overlapping genes associated with gray matter atrophy in PD patients are illustrated in Fig. 3. These genes were exclusively expressed in the brain tissue by tissue-specific expression analysis (Fig. 3a). Cell type-specific expression analysis demonstrated that these genes were expressed specifically among dopamine receptor cells (caudate putamen D1 and D2 dopamine receptor) and multiple types of neurons, including the cortex, corticothalamic neurons, caudate putamen, corticopontine neurons, corticosterone-expressing neurons, and prepronociceptin-expressing neurons (Fig. 3b). Regarding the temporal specific expression, it was discovered that these genes were expressed throughout almost the entire developmental stage (Fig. 3c).
Specific expression analyses of the identified genes associated with gray matter atrophy in PD patients. a) Tissue-specific expression. The radial array numbers are −log10(q). b) Cell-specific expression. The x axis indicates −log10(q), and the y axis indicates cell types. c) Temporal specific expression. The x axis denotes developmental stages and the y axis denotes −log10(q). The gray dashed line represents the statistical significance threshold of q < 0.05. Abbreviations: Cort+, corticosterone-expressing neurons; Cpu, caudate putamen; Cpu.D1, caudate putamen D1 dopamine receptor; Cpu.D2, caudate putamen D2 dopamine receptor; Ctx, cortex; Glt25d2, corticopontine neurons; Ntsr+, corticothalamic neurons; Pnoc+, prepronociceptin-expressing neurons; pSI, specificity index probability.
PPI network and hub genes
According to PPI analysis, an interconnected PPI network could be created by 742 genes from the 1,952 genes (Fig. 4a). This network consisted of 919 edges, which was significantly higher than expected (P = 1.47 × 10−8). Sixteen genes with the top 10% highest degree values were defined as hub genes (Supplementary Table 5). In addition, we delineated the spatiotemporal expression trajectory of three hub genes with the highest degree values (e.g., CTNNB1, MAPK3, and CALM3) (Fig. 4b).
The protein-protein interaction (PPI) network and hub genes. a) PPI network with 742 genes and 919 edges constructed by the identified genes associated with gray matter atrophy in PD. Three most contributive hub genes (e.g., CTNNB1, MAPK3, and CALM3) are labeled by boxes. The P-value denotes the statistical significance of how likely the proteins could construct a network. b) Spatial temporal expression curves of three representative hub genes (e.g., CTNNB1, MAPK3, and CALM3) with the highest degree values.
Behavioral relevance analysis
We discovered that these genes associated with gray matter atrophy in PD were correlated with an array of behavioral terms, including action (e.g., execution, motor learning, and inhibition), cognition (e.g., spatial and somatic cognition, implicit memory, language orthography), intensity emotion, gastrointestinal-genitourinary, and perception (e.g., motion, shape, and color vision perception) by correlating gene expression with behavioral domains using BrainMap (Fig. 5a). In the same way, we found that the identified genes were related to behavioral domains including visual and early visual, motion, emotion, attention, and dementia (Fig. 5b) via Neurosynth, which were consistent with those from the BrainMap.
The behavioral relevance of the identified genes associated with gray matter atrophy in PD patients via BrainMap a) and Neurosynth b). The size and frequency of the behavioral domain terms in the word cloud diagrams represent the correlation between expression measures of the identified genes and activation values. Abbreviations: GI-G, gastrointestinal-genitourinary.
Discussion
In summary, we explored a new approach to shed light on the genetic mechanisms underlying gray matter atrophy in PD by conducting a combined analysis of brain imaging and gene expression data. The neuroimaging meta-analysis, along with the independent morphology analysis, altogether showed that PD patients consistently exhibit significant gray matter atrophy in the superior temporal gyrus. Furthermore, a spatial correlation study between the transcriptome data and neuroimaging phenotype indicated that the identified gray matter atrophy was spatially related to the expression of 1,952 overlap genes, which were enriched for a diverse range of MFs, BPs, and CCs as well as various biological pathways. In addition, these genes were found to be selectively expressed in brain tissue, dopamine receptor cells, and neurons throughout almost the entire developmental stage. Likewise, these genes showed the potential for creating a PPI network supported by 16 putative hub genes of functional significance. Overall, our findings indicate that gray matter atrophy in PD could potentially be a consequence of intricate interactions among a complex set of genes, confirming the polygenic nature of this neurological condition.
To characterize GMV changes while balancing generalizability and data authenticity, we conducted a neuroimaging meta-analysis alongside a morphological study using the VBM method in an independent dataset. This approach not only advances our understanding of PD-related phenotypes but also ensures that the subsequent transcriptome-neuroimaging association analysis is not entirely influenced by the indirect effect values. This pipeline has been adopted in previous studies, demonstrating its stability and reliability (Sun et al. 2023; Cai et al. 2024). Specifically, the updated neuroimaging meta-analysis indicated significant gray matter atrophy in the bilateral superior temporal gyrus and inferior parietal gyrus, the right angular gyrus, and the left superior frontal gyrus in PD patients. Meanwhile, PD patients from the independent dataset also exhibited gray matter atrophy in the left superior temporal gyrus and amygdala, plausibly consistent with findings from the meta-analysis. The overlap results emphasize the significant appearance of gray matter atrophy in the superior temporal gyrus, consistent with previous research findings (Xu et al. 2020; Zheng et al. 2022; Zhong et al. 2022; Baagil et al. 2023). Abnormalities in this region have been linked with key symptoms of PD, including theory of mind impairment, apathy, cognitive impairment, depression, freezing of the gait, and frequent falls (Lee et al. 2013; Otomune et al. 2019; Li et al. 2020; Orso et al. 2020). Furthermore, previous research has confirmed that Lewy bodies and Lewy neurites can also be detected outside the substantia nigra in PD (Braak et al. 2003). Braak and colleagues have proposed the classic six stages of PD progression, starting in the peripheral nervous system and progressively affecting the central nervous system in a caudal-to-rostral direction within the brain (Braak et al. 2003). With the progression of PD, alpha-synuclein inclusions collect in the temporal cortex with some preference, and PD patients lose dopamine receptors in the temporal lobe (Goedert et al. 2013; Picco et al. 2015). Understanding the patterns of GMV loss may provide valuable insights into the pathogenesis associated with aging and neurodegeneration, which could be useful for developing future interventions and studying cognitive reserve and compensation. However, the inconsistencies should not be ignored, as the right hemisphere exhibited much more apparent atrophy compared to the left in the meta-analysis, whereas the independent VBM analysis showed observable results only in the left hemisphere. Though the transcription–neuroimaging association analysis was minimally affected by the inconsistencies in the right hemisphere, since tissue samples only from the left cerebral cortex were used in this study—all six donors had expression data for the left hemisphere, while only two had samples for the right (Hawrylycz et al. 2012). Some disparities could be interpreted from various perspectives, such as the relatively limited sample size, milder severity of the disease, and shorter disease duration of the PD patients in our independent dataset compared to previous studies. These limitations highlight the need to expand our participant dataset in future research.
In particular, we discovered 1,952 genes whose expression patterns were associated with gray matter atrophy by means of transcription-neuroimaging spatial association analysis. When compared with the 90 independent genome-wide significant association signals identified from the large meta-analysis, we found 14 overlap genes, including the classic disease-causing genes (SNCA and LRRK2) of PD (Deng et al. 2018). Our surprising discovery of gene expression can be linked to the detoxification of heavy metals. It is generally accepted that the core motor impairments in PD are triggered by the degeneration of dopaminergic neurons in the substantia nigra. These neurons are autonomous pacemakers with large cytosolic Ca2+ oscillations that have been related to basal mitochondrial oxidant stress and turnover (Zampese and Surmeier 2020). More importantly, potassium channels are the most prevalent ion channel family and have been proven to be of importance in PD pathogenesis due to their functions in influencing neuronal excitability, neurotransmitter release, synaptic transmission, and neuroinflammation (Chen et al. 2023). In terms of the role of disturbed synaptic function in PD, numerous genes altered in familial types of PD are related to systems of synaptic homeostasis and synaptic function (Kumaran and Cookson 2015), and genes implicated in sporadic PD also contribute to these processes (Soukup et al. 2018). The spatial associations between gene expression and neuroimaging phenotypes of diseases are under active investigation, which may help to further elucidate the genetic mechanisms underlying brain abnormalities that characterize diseases.
Specific expression analyses revealed that the genes associated with gray matter atrophy in PD were specifically expressed in brain tissue, dopamine receptor cells, and various neurons. The D1- and D2- dopamine receptors, which signal through Gs and Gi, respectively, represent the principal stimulatory and inhibitory dopamine receptors in the central nervous system (Zhuang et al. 2021). The dopamine D2 receptor, especially found on CD4+ T cells, is typically thought to be protective against neuroinflammation and neurodegeneration in PD (Liu et al. 2021). It has also become apparent that an ideal antiparkinsonian medication would need to stimulate both the D1 and D2 dopamine receptors, which represent the primary therapeutic targets (De Keyser et al. 1995). It is also worth noting that these identified genes have the potential to create a PPI network supported by 16 putative hub genes with functional value in comprehending the pathology and therapy of PD. CTNNB1, MAPK3, and CALM3 contributed most significantly to the PPI network as hub genes. The reduction in protein expression of the neurogenesis-related gene CTNNB1 in midbrain dopaminergic neurons demonstrated a significant relationship between the epigenetic dysregulation of Wnt signaling and the pathophysiology and development of PD (Zhang et al. 2016). In vivo experiments confirmed that Nootkatone could inhibit the expression of MAPK3 by activating the PI3K/Akt signaling pathway, reducing neuroinflammation, and ultimately ameliorating rotenone-induced PD symptoms (Yao et al. 2022). MAPK3 has also been prioritized as a potential target for neurodegenerative diseases and has provided preliminary evidence for drug development (Ge et al. 2023). Additionally, CALM3 might also play an important role in PD pathogenesis (George et al. 2019). Regarding behavioral relevance, beyond the universally recognized movement disorders, cognitive impairment, and dementia, visual symptoms such as trouble reading, double vision, illusions, and complex visual hallucinations are also common in PD and dementia (Archibald et al. 2013). Moreover, gastrointestinal and genitourinary diseases may appear decades before PD diagnosis, highlighting their potential as sites for establishing early diagnostic testing and studying the disease’s pathophysiology (Scott et al. 2021). Traditionally, PD was believed to primarily impact the motor system due to the loss of dopaminergic neurons in the substantia nigra. However, recent studies have revealed that PD could be attributed to abnormal aggregation of α-synuclein and the spreading of pathology between the gut, brainstem, and higher brain regions underlying the disease’s development and progression (Morris et al. 2024). Specifically, PD appears to result from the complex interplay of aberrant α-synuclein aggregation, mitochondrial dysfunction, lysosomal or vesicle transport issues, synaptic transport problems, and neuroinflammation, which collectively result in accelerated neuronal death of primarily dopaminergic neurons (Kalia and Lang 2015).
Several limitations should be considered when interpreting our findings. First, we merely used left hemispheric data as gene expression data as the right hemisphere was only accessible from two donor brains (Hawrylycz et al. 2012; Arnatkeviciute et al. 2019). However, the results from the meta-analysis showed more reductions in the right hemisphere compared with the left. Second, the gene expression data came from post-mortem brains, whereas the neuroimaging data came from living brains. To address potential variability, we set a DS threshold to focus on genes with more consistent expression patterns, making it feasible to conduct the transcription–neuroimaging spatial association analysis (Zeng et al. 2012). Third, since only the peak coordinates and related effect sizes of significant clusters reported in prior research were employed, the meta-analysis z map cannot capture the nature and extent of whole-brain gray matter differences between PD patients and HC. Fourth, while the neuroimaging results from the independent dataset align with parts of the meta-analysis, the balance of sample size between them was worthy of further promotion. Despite setting strict corrections (FWE correction and standard random-effects model) to increase the rigor of the meta-analysis, we acknowledge the need to expand the size and variety of our independent dataset in future studies. Lastly, our study is a general analysis of PD without subdivision into specific subtypes or clinical phenotypes. Studies focusing on the progression and clinical relevance of brain gray matter atrophy in PD patients will be warranted in the future.
In conclusion, our data revealed prominent gray matter atrophy in PD by combining a comprehensive meta-analysis and an independent dataset VBM analysis. Additionally, we discovered that the identified gray matter atrophy was spatially associated with the expression levels of 1,952 genes characterized by a variety of functional characteristics. Our findings may not only offer unique insight into the genetic mechanisms of gray matter atrophy in PD but also inform potential treatment approaches targeting the molecular substrates underlying the brain morphological abnormalities of this disorder.
Author contributions
Yi Ji (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing—original draft), Min Xu (Conceptualization, Data curation, Resources, Validation), Han Zhao (Conceptualization, Data curation, Formal analysis, Resources, Software, Supervision, Validation), Huanhuan Cai (Conceptualization, Data curation, Formal analysis, Resources, Software, Supervision, Validation), Kaidong Chen (Formal analysis, Investigation, Validation, Visualization), Li Zhang (Data curation, Resources), Haixia Mao (Investigation, Supervision), Feng Wang (Data curation, Investigation, Resources, Supervision), Jiajia Zhu (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing—review & editing), and Xiangming Fang (Conceptualization, Funding acquisition, Project administration, Supervision, Writing—review & editing). All authors have read and approved the final manuscript. Y.J. and M.X. share the first authorship of this article.
Funding
This work was supported by Medical Expert Team Program of Wuxi Taihu Talent Plan (THRC-TD-YXYXK-2021), Wuxi Medical Innovation Team Program (CXTD2021002), Natural Science Foundation of Jiangsu Province (No. BK20191143, X.F.), and National Natural Science Foundation of China (No. 81271629, X.F.).
Conflict of interest statement: All authors declare no financial or non-financial competing interests to disclose regarding this article.
Data availability
The brain gene expression data are publicly available from the open Allen Human Brain Atlas (AHBA) dataset (http://www.brain-map.org). The individual de-identified participant data from the independent dataset will be shared by the corresponding author upon reasonable request.
Ethics approval and patient consent statement
The study received official approval (KY21133) from the Ethics Committee of the Affiliated Wuxi People’s Hospital of Nanjing Medical University, and all participants provided written informed consent before participation.
References
Author notes
Yi Ji and Min Xu share the first authorship of this article.
