Genetic risk for Alzheimer’s dementia predicts motor deficits through multi-omic systems in older adults

Alzheimer’s disease manifests with both cognitive and motor deficits. However, the degree to which genetic risk of Alzheimer’s dementia contributes to late-life motor impairment, and the specific molecular systems underlying these associations, are uncertain. Here, we adopted an integrative multi-omic approach to assess genetic influence on motor impairment in older adults and identified key molecular pathways that may mediate this risk. We built a polygenic risk score for clinical diagnosis of Alzheimer’s dementia (AD-PRS) and examined its relationship to several motor phenotypes in 1885 older individuals from two longitudinal aging cohorts. We found that AD-PRS was associated with a previously validated composite motor scores and their components. The major genetic risk factor for sporadic Alzheimer’s dementia, the APOE/TOMM40 locus, was not a major driver of these associations. To identify specific molecular features that potentially medicate the genetic risk into motor dysfunction, we examined brain multi-omics, including transcriptome, DNA methylation, histone acetylation (H3K9AC), and targeted proteomics, as well as diverse neuropathologies. We found that a small number of factors account for the majority of the influence of AD-PRS on motor function, which comprises paired helical filament tau-tangle density, H3K9AC in specific chromosomal regions encoding genes involved in neuromuscular process. These multi-omic factors have the potential to elucidate key molecular mechanisms developing motor impairment in the context of Alzheimer’s dementia.

assessed in the substantia nigra in the mid to rostral midbrain, near or at the exit of the third cranial nerve, using H&E stain and six micron sections using a dichotomous scale (absent or present) 6 . Lewy bodies were assessed in six regions using a monoclonal phosphorylated antibody to a-synuclein 7 . TDP-43 staging from amygdala to limbic and neocortical regions was determined using monoclonal TDP-43 antibody 8 . Hippocampal sclerosis was evaluated in a coronal section of the mid-hippocampus at the level of the lateral geniculate body and was reported as present if there was severe neuronal loss and gliosis in the CA1 sector, subiculum or both 9 . Chronic macroscopic and microinfarcts infarcts were recorded during gross examination and confirmed histologically 10 . Cerebral amyloid angiopathy (CAA) was assessed in four neocortical regions using immunohistochemistry 11 . β-amyloid depositions in meningeal and parenchymal vessels in each region were rated and averaged to obtain a summary CAA measure.
Severity of atherosclerosis was graded by gross examination of vessels in the circle of Willis, and arteriolosclerosis was graded on H&E stained sections of basal ganglia 12 . The complete list of brain pathologies assessed in this study is in Supplementary Table 2.

Omics measurements
Details on omics data processing were published previously 13 . Briefly, genotyping data was measured using DNA extracted from peripheral blood mononuclear cells or frozen brain tissue, and quality control steps were performed as described previously 13,14 . After the quality control steps, genotyping data included 7,159,943 SNPs in 2,093 subjects. DNA methylation data were generated using DNA extracted from DLPFC 15 and pre-processed as described previously 16 . The pre-processed data consisted of ~130,000 methylation loci in 533 subjects. Histone H3 acetylation on lysine 9 (H3K9AC) data was measured by Chromatin Immunoprecipitation (ChIP) assay using anti-H3K9AC mAb coupled with sequencing was performed in gray matter of DLPFC 17 . Quality control steps were described previously 13,17,18 . The pre-processed data consisted of 26,384 histone peaks in 516 subjects. RNAseq data was generated from DLPFC and quality control steps were performed as described previously 19,20 . The pre-processed data consisted of 13,484 genes in 432 subjects. To alleviate a large multiple testing burden, we followed the standard practice of reducing DNA methylation, histone acetylation, and gene expression to comethylated, coacetylated, and coexpressed modules 21 , each of which was composed of variables with similar patterns of methylation, acetylation, or expression, as measured across all individuals. Using SpeakEasy 22 , we identified 58 DNA comethylation modules, 80 histone coacetylation modules, and 49 coexpression modules 13 . The miRNA expression profiles were measured in DLPFC samples using the NanoString nCounter miRNA expression assay and pre-processing steps were performed as described previously 13,23 . The preprocessed data consisted of 292 miRNAs in 543 subjects. Selected reaction monitoring (SRM) proteomics was performed using frozen DLPFC tissue for 67 proteins selected by the consortium members of Accelerating Medicines Partnership for Alzheimer's Disease (AMP-AD; https://www.synapse.org/#!Synapse:syn2580853) 24 . Pre-processing steps were described previously 13,24 . The pre-processed data consisted of 67 proteins in 527 subjects.

Gene ontology (GO) enrichment analysis for histone coacetylation modules
To examine whether histone coacetylation modules were involved in particular biological processes, we performed GO enrichment analysis. GO gene sets were downloaded from MSigDB v6.1 25,26 . The GREAT algorithm 27

Estimation of AD-PRS effect explained by endophenotypes
To evaluate the proportion of AD-PRS effect on motor function explained by endophenotypes, we compared the variance of motor function explained by AD-PRS and that given each molecular phenotype as follows: First, the total variance of motor function was computed as the total sum of squares (SS), where . is an intercept, . is a coefficient for AD-PRS, . is a coefficient for a mediator, and .
is an error term. The SS for the residual of the second model was computed as 1 * 100.

Estimation of Bayesian network
To infer the relationships among AD-PRS, endophenotypes, and a motor function, we used a Bayesian network, which is a multivariate probabilistic model whose conditional independence relations can be represented by a directed acyclic graph (DAG) with vertices = P 1 , … , R S, and directed edges ( , ) ∈ ⊂ × (note that we use the notation and * , interchangeably, to refer to a node). A vertex in a DAG corresponds to a random variable Z in the Bayesian network. Assuming the local directed Markov property, each variable is independent of its nondescendant variables conditional on its parent variables. Thus, the state of Z can be determined only by the state of parent variables, which is formally expressed by the conditional probability, where Z state occurs under given parents' state ]^. Therefore, the probability where observed data, , is generated from a given DAG can be factored as ( | ) = , where = P 1 , … , R S a , Z is the set of parents of , and ]^ = b * : ∈ Z d.
To learn the DAG structure, which is the process of finding with high ( | ), we used a where ( ) is a prior on the network structure , and represents the space of all DAGs with vertices. The MCMC sampling allows us to obtain ensembles of DAGs with high ( | ) and avoid overfitting to the data. To utilize genetic information as a clue to infer the directions of other edges, we restricted a direction of edges so that AD-PRS can have only out-going edges to other nodes. We did not pose any restrictions for non-genetic nodes. We ran 300,000 steps of MCMC sampling using the REV algorithm 31 and discarded the first 2,400 steps as a burn-in.
Then, edge frequencies in the sampled networks were counted, and a consensus network was generated by taking the regulation that presented most frequently among the three possible: node1 regulates node2, node2 regulates node1, and node1 is independent of node2. The detailed implementation of learning network structure based on systems genetics data is described elsewhere 32 . Given the estimated network structure, each variable was regressed with its parent variables and obtained variance explained by each parent variable and p-value associated with it.