阿尔茨海默病脑和脑脊液的大规模蛋白质组学分析揭示了与小胶质细胞和星形胶质细胞活化相关的能量代谢的早期变化

共有AD蛋白共表达网络的构建与验证

We analyzed dorsolateral prefrontal cortex (DLPFC) tissue in 44 cases from the Baltimore Longitudinal Study of Aging (BLSA), 178 cases from the Banner Sun Health Research Institute (Banner), 166 cases from the Mount Sinai School of Medicine Brain Bank (MSSB), and 65 cases from the Adult Changes in Thought Study (ACT), for a total of 453 control, asymptomatic AD (AsymAD), and AD brains (Figure 1A, Supplementary Table 1). AsymAD was defined as postmortem pathology consistent with an AD diagnosis but without dementia, based on the NIA research framework for AD². Tissues were analyzed by mass spectrometry-based proteomics using label-free quantitation (LFQ), and the resulting mass spectrometry data were processed using a common pipeline to arrive at 5688 total quantified proteins. We included proteins with fewer than 50% missing values in the subsequent analyses, as it was determined that this threshold was robust to potential spurious correlations given the power of the study (Extended Data Figure 1). We also removed by regression the effects of age, sex, and post-mortem interval on the protein quantitative data, even though these covariates did not strongly influence the data (Extended Data Figure 2). The final adjusted 3334 proteins were used to generate a protein co-expression network using the weighted correlation network analysis (WGCNA) algorithm. The resulting network consisted of 13 protein co-expression “modules,” or communities of proteins with similar expression patterns across the cases analyzed (Figure 1B, Supplementary Figures 1 and 2, Supplementary Table 2). The modules ranged in size from 254 proteins (M1) to 20 proteins (M13). These modules could also be identified independently of the WGCNA algorithm using t-distributed stochastic neighbor embedding (t-SNE) analysis (Extended Data Figure 3), demonstrating that the protein communities identified by the WGCNA algorithm were robust. Gene ontology (GO) analysis of the protein module members revealed a clear ontology for eleven out of the thirteen modules, encompassing a diverse mix of biological functions, processes, and components (Figure 1B, Supplementary Figure 3). To assess whether a given co-expression module was related to AD, we correlated the module eigenprotein—or first principle component of the module protein expression level—to the neuropathological hallmarks of AD: amyloid-β plaques and neurofibrillary tangles. We also correlated the module eigenproteins to cognitive function as assessed by the Mini-Mental Status Examination (MMSE), and functional status as assessed by the Clinical Dementia Rating Scale (CDR), at the last research evaluations prior to death to capture module-disease relationships that may be independent of amyloid-β plaque or tau tangle pathology (Figure 1B, Supplementary Figure 1). We observed six modules that were significantly correlated with all pathological, cognitive, and functional measures, and whose ontologies could be best characterized by a structural component or a biologic process: modules M1 synapse, M3 mitochondrial, M4 glucose and carbohydrate metabolism (subsequently referred to as sugar metabolism), M5 extracellular matrix, M6 cytoskeleton, and M10 RNA binding/splicing. The M4 sugar metabolism module showed the strongest AD trait correlations (cognition r=−0.67, p=8.5e⁻²³; neurofibrillary tangle r=0.49, p=4.7e⁻²⁷; amyloid-β plaque r=0.46, p=1.3e⁻²³; functional status r=0.52, p=2.6e⁻¹²). Because AD neuropathology is not homogenous even within the same brain region, and because neuropathological measurements of AD pathology are semi-quantitative and subject to a certain degree of individual variability in assessment¹⁵, we also correlated module eigenproteins to mass spectrometry measurements of amyloid-β and the tau microtubule binding region, which comprises neurofibrillary tangles, within the DLPFC tissue used for proteomic analysis (Extended Data Figure 4). We observed strong concordance between neuropathological and molecular measurements of AD pathology. 重试    错误原因

Because many protein co-expression changes in the brain can be driven by cell type changes^16,17, we also assessed the cell type nature of each co-expression module by asking whether the module was enriched in particular cell type marker proteins (Figure 1B). We observed significant enrichment of neuronal proteins in the M1 synapse module and enrichment of oligodendrocyte markers in the M2 myelin module, as expected. We also observed enrichment of astrocyte and microglial proteins in the M4 sugar metabolism module, microglial and endothelial proteins in the M5 extracellular matrix module, and endothelial markers in the M7 translation/ribosome module. These findings suggest that the biological processes reflected by GO analysis for each module may be altered in AD within a particular cell type. To incorporate the cell type nature of each module into its description, we subsequently refer to those modules with strong cell type enrichment as the “M1 synapse/neuron” module, the “M2 myelin/oligodendrocyte” module, the “M4 astrocyte/microglial metabolism” module, and the “M5 endo/micro extracellular matrix” module.

To assess the relationship of the network modules to diagnostic classification, we measured the module eigenprotein values by case status (Figure 1C and Supplementary Figure 1). In general, most modules that were increased or decreased in AD compared to control also showed a trend, or were significantly changed, in the same direction in the AsymAD group, indicating that these modules reflect pathophysiologic processes that begin early—in the preclinical phase—of AD. The M1 synapse/neuron, M3 mitochondrial, and M4 astrocyte/microglial metabolism modules showed the strongest differences by case status. We also assessed the influence of APOE genotype on module eigenproteins, but did not find strong effects except for the APOE ε2 allele on the M2 myelin/oligodendrocyte module, which appeared to attenuate the observed changes in AD (Supplementary Figure 1).

To validate the AD network, we analyzed 340 DLPFC brain tissues from a community-based aging cohort, the Religious Orders Study and Memory and Aging Project (ROS/MAP)^18–20, with a different mass spectrometry-based protein quantification approach using isobaric multiplex tandem mass tags (TMT)^21–23. A protein co-expression network was constructed from the ROS/MAP cases, and network module preservation statistics as well as synthetic module eigenproteins were used to assess conservation of the consensus AD LFQ-based network in the ROS/MAP TMT-based network (Extended Data 5B and C, Supplementary Figure 4). We found that all consensus LFQ modules were preserved in the ROS/MAP TMT-based network. Furthermore, targeted protein measurements in a cohort of 1016 ROS/MAP control, AsymAD, and AD brains by another mass spectrometry protein quantification approach—selected reaction monitoring (SRM)—showed that individual module proteins had the same direction of change as the AD LFQ-based network co-expression module of which they were a member (Supplementary Figures 5 and 6). In summary, we were able to construct a robust AD protein co-expression network from mass spectrometry-based proteomic analysis of greater than 450 human DLPFC brain tissues from multiple centers. We found that many of these modules correlated with AD neuropathology and cognitive function, reflected a number of different biological processes and cell types, and were altered in the preclinical stage of AD.

The AD Network Is Preserved in Other Brain Regions

The consensus AD network was generated from analysis of DLPFC tissue. To assess whether the network was similar in other brain regions commonly affected in AD, we analyzed control and AD brain tissue from temporal cortex in a separate set of 111 brains from the Mayo Clinic, and control, AsymAD, and AD brain tissue from precuneus in the same set of brains from the BLSA (Figure 2A) using LFQ-MS. Co-expression networks were built for each brain region, and network preservation statistics were used to assess module preservation from DLPFC in temporal cortex (Figure 2B) and precuneus (Figure 2C). We found that all consensus AD network modules derived from DLPFC were preserved in temporal cortex, and twelve out of the thirteen modules were preserved in precuneus. Analysis of synthetic module eigenprotein values by case status showed similar differences between and among case groups in temporal cortex (Figure 2D, Supplementary Figure 7) and precuneus (Figure 2E, Supplementary Figure 8) brain regions, with changes in AsymAD more pronounced in precuneus than in DLPFC. These findings suggest that the consensus AD network is generalized across brain regions that are commonly affected in AD.

Effects of Aging on AD Network Modules

To better understand the influence aging—the strongest risk factor for AD—may have on the consensus AD network, we analyzed DLPFC tissues from Johns Hopkins in 84 cases ages 30 to 69 (Figure 3A) by LFQ-MS. All cases had a final primary neuropathological diagnosis of control. We created synthetic eigenproteins in the aging cohort from the consensus AD network modules and asked whether the synthetic module eigenproteins changed with age (Figure 3B and Supplementary Figure 9). We found that the M1 synapse/neuron and M4 astrocyte/microglial metabolism modules decreased and increased with aging, respectively, while the M3 mitochondrial and M10 RNA binding/splicing modules were not affected by aging. Other modules that appeared to be affected by aging included the M6 cytoskeleton, M7 translation/ribosome, and M9 translation/ribosome modules (Supplementary Figure 9). Additional information on the correlation of individual proteins with age and overlap with markers of cellular senescence is provided in Supplementary Table 3. These findings indicate that the relationship between aging and AD at the proteomic level is complex, and that some, but not all, AD trait-associated modules are influenced by the aging process.

AD Network Changes in Other Neurodegenerative Diseases

To explore the specificity of these network changes for AD, we analyzed 331 DLPFC tissues by LFQ-MS from control, AD, amyotrophic lateral sclerosis (ALS), frontotemporal lobar degeneration with TDP-43 pathology (FTLD-TDP), progressive supranuclear palsy (PSP), corticobasal degeneration (CBD), Parkinson’s disease and Parkinson’s disease dementia (PD/PDD), and multiple systems atrophy (MSA) cases (Figure 4A). We created synthetic eigenproteins for consensus AD network modules and assessed whether they changed in these different neurodegenerative diseases compared to AD (Figure 4B, Supplementary Figure 10, Supplementary Table 4). We found that the M1 synapse/neuron and M4 astrocyte/microglial metabolism modules showed significant changes in FTLD-TDP and CBD cases, similar to AD, whereas the M3 mitochondrial and M10 RNA binding/splicing modules showed more mixed changes across other diseases. To further validate these findings, we used a targeted mass spectrometry method called parallel reaction monitoring (PRM)²⁴ to measure 323 individual proteins from approximately one-third of the cases analyzed in the untargeted experiments (Supplementary Figure 11, Supplementary Table 4). Protein levels across all cases were highly correlated (r=0.92, p=1.4e⁻¹²⁴) between LFQ and PRM measurements (Extended Data Figure 6B). We created synthetic eigenproteins from these targeted PRM protein measurements by AD consensus module, and assessed eigenprotein changes by disease category (Extended Data Figure 6C, Supplementary Figure 12, Supplementary Table 4). We observed very similar AD network module changes across diseases compared to the untargeted measurements, validating the findings from the untargeted LFQ measurements. These results indicate that certain AD network modules are affected to a greater extent in AD compared to other neurodegenerative diseases, and that FTLD and CBD show many similar changes to AD, with the caveat that not all neurodegenerative diseases affect the DLPFC region equally at end-stages of disease.

The M4 Astrocyte/Microglial Metabolism Module is Enriched in AD Genetic Risk Factors and Markers of Anti-Inflammatory Disease-Associated Microglia

We applied an algorithm to calculate a weighted disease risk score for proteins according to their linkage disequilibrium with AD-associated single nucleotide polymorphisms (SNPs) discovered through AD genome wide association studies (GWAS)²⁵. We then calculated whether a given AD network module was enriched in these risk factor proteins. We found that the M2 myelin/oligodendrocyte and M4 astrocyte/microglial metabolism modules were significantly enriched in gene products contained within AD risk factor loci (Figure 5A), suggesting that the biological functions or processes reflected by these protein co-expression modules may serve causative roles in AD.

Given the strong AD trait associations of the M4 astrocyte/microglial metabolism module and its enrichment in AD genetic risk factors, we more deeply investigated the cell type nature of this co-expression module. Although expression of the M4 astrocyte/microglia metabolism module is increased with progression from a normal to an AD disease state, and a majority of the most significantly increased proteins in AD are members of this module (Extended Data Figure 7), it is unclear whether these glial responses are deleterious or protective. To better understand the role of these glial cell type responses in AD, we first examined differential expression of astrocyte and microglia protein markers in AD brain by the types of cellular phenotypes with which they are associated in AD animal models^26–31. We found that for both astrocytic markers (Extended Data Figure 8) and microglial markers (Extended Data Figure 9), there appeared to be a bias towards expression of markers that are generally considered to be protective. We formally tested this observation with marker over-representation analysis in the AD network (Figure 5B, Supplementary Table 5). Microglial protein markers that are increased in response to amyloid-β plaques but decreased in response to LPS—or markers of anti-inflammatory disease-associated microglia²⁷—were significantly enriched in the M4 module. Astrocyte markers were more mixed in module M4, with a majority of markers being shared between deleterious A1 and protective A2 phenotypes²⁶. Astrocyte and microglia phenotype markers that overlap with the top 100 proteins by module eigenprotein correlation value in the M4 module are shown in Figure 5C. The majority of these markers were from microglia (Supplementary Table 5). To further validate these findings, we analyzed whether these markers were increased at both the transcript and protein levels in acutely isolated microglia from AD mouse models^32,33. The top 30 most differentially abundant microglial transcripts corresponding to proteins in the M4 module were found to be heavily biased toward an anti-inflammatory phenotype (Figure 5D, Supplementary Table 5). Furthermore, many of the disease-associated M4 microglial protein markers were found to be increased in microglia undergoing active amyloid plaque phagocytosis (Extended Data Figure 10, Supplementary Table 5)³³. In summary, we found that the M4 astrocyte/microglial metabolism module was enriched in AD genetic risk factors, and that microglia cell type markers within M4 appeared to be biased towards a protective anti-inflammatory, rather than a deleterious pro-inflammatory, microglial phenotype.

M4 Astrocyte/Microglial Metabolism Proteins Are Increased in Cerebrospinal Fluid

To explore whether proteins from the M4 astrocyte/microglial metabolism module might also be able to serve as AD fluid biomarkers, we analyzed CSF from two separate cohorts: one cohort of 297 subjects consisting of controls and AD patients (Cohort 1), and a second cohort of 96 subjects classified into control, AsymAD, and AD (Cohort 2). Subjects in both cohorts were classified by the “A/T/N” AD biomarker classification framework (Figure 6A)³⁴. CSF from both cohorts was analyzed using a TMT-MS approach without prior pre-fractionation and without depletion of highly abundant proteins. In Cohort 1, we observed 22 proteins that mapped to the M4 astrocyte/microglial metabolism module in brain (Supplementary Figure 13). All of them showed either an increase in AD or no change, with 10 reaching statistical significance at p < 0.05. The most significantly increased M4 module proteins observed in Cohort 1 are shown in Figure 6B, and include the M4 hub proteins CD44, peroxiredoxin-1 (PRDX1), and dimethylarginine dimethylaminohydrolase-2 (DDAH2), in addition to the metabolic proteins lactate dehydrogenase B-chain (LDHB) and pyruvate kinase (PKM) involved in glycolysis. To validate these findings, and to assess whether the observed changes in CSF levels of M4 proteins occur prior to the development of cognitive impairment, we analyzed subjects in Cohort 2, approximately one-third of which had AsymAD. AsymAD was defined as CSF levels of amyloid-β, total tau, and phospho-tau consistent with an AD diagnosis, but without cognitive impairment. In Cohort 2, 27 proteins mapped to the M4 astrocyte/microglial metabolism module in brain (Supplementary Figure 14). Of these 27, 17 overlapped with M4 proteins measured in discovery Cohort 1, and showed the same direction of change in AD CSF. In addition, many also showed significant or trend elevations in AsymAD, including CD44, LDHB, and PKM, and correlated with cognitive function (Figure 6C). In summary, multiple M4 astrocyte/microglial metabolism module protein members could be measured in human CSF by mass spectrometry without fractionation or prior depletion of highly abundant proteins. A number of these proteins were elevated in AsymAD and AD, including M4 hub proteins CD44, PRDX1, and DDAH2.

Brain Tissue Samples and Case Classification

Brain tissue used in this study was obtained from the autopsy collections of the Baltimore Longitudinal Study of Aging⁵⁵, Banner Sun Health Research Institute⁵⁶, Mount Sinai School of Medicine Brain Bank, Adult Changes in Thought Study, Mayo Clinic Brain Bank, Religious Orders Study and Rush Memory and Aging Project⁵⁷, University of Pennsylvania School of Medicine Brain Bank, and the Baltimore Coroner’s Office. Tissue was from the dorsolateral prefrontal cortex (Brodmann Area 9 where available), or temporal cortex and precuneus regions where indicated. Human postmortem tissues were acquired under proper Institutional Review Board (IRB) protocols at each respective institution. Postmortem neuropathological evaluation of neuritic plaque distribution was performed according to the Consortium to Establish a Registry for Alzheimer’s Disease (CERAD) criteria¹⁵, while extent of spread of neurofibrillary tangle pathology was assessed with the Braak staging system⁵⁸. Other neuropathologic diagnoses were made in accordance with established criteria and guidelines⁵⁹. All case metadata, including age, sex, post-mortem interval, cognitive function, APOE genotype, neuropathological criteria, and disease status, are provided in Supplementary Table 1. Case classification harmonization across cohorts was performed using the following rubric: cases with CERAD 0-1 and Braak 0-3 without dementia at last evaluation were defined as control (if Braak equals 3, then CERAD must equal 0); cases with CERAD 1-3 and Braak 3-6 without dementia at last evaluation were defined as AsymAD; cases with CERAD 2-3 and Braak 3-6 with dementia at last evaluation were defined as AD. Dementia was defined as MMSE <24, CASI score <81, or CDR ≥1, based on prior comparative study⁶⁰. Mayo and UPenn cases were not included in the case harmonization scheme, and therefore preservation of consensus network modules in these cohorts provides an additional degree of robustness.

Brain Tissue Homogenization and Protein Digestion

Procedures for tissue homogenization for all tissues were performed essentially as described^17,21. Approximately 100 mg (wet tissue weight) of brain tissue was homogenize in 8 M urea lysis buffer (8 M urea, 10 mM Tris, 100 mM NaHPO4, pH 8.5) with HALT protease and phosphatase inhibitor cocktail (ThermoFisher) using a Bullet Blender (NextAdvance). Each Rino sample tube (NextAdvance) was supplemented with ~100 μL of stainless steel beads (0.9 to 2.0 mm blend, NextAdvance) and 500 μL of lysis buffer. Tissues were added immediately after excision and samples were then placed into the bullet blender at 4 °C. The samples were homogenized for 2 full 5 min cycles, and the lysates transferred to new Eppendorf Lobind tubes. Each sample was then sonicated for 3 cycles consisting of 5 s of active sonication at 30% amplitude, followed by 15 s on ice. Samples were then centrifuged for 5 min at 15,000 x g and the supernatant transferred to a new tube. Protein concentration was determined by bicinchoninic acid (BCA) assay (Pierce). For protein digestion, 100 μg of each sample was aliquoted and volumes normalized with additional lysis buffer. For the ROS/MAP cohort, an equal amount of protein from each sample was aliquoted and digested in parallel to serve as the global pooled internal standard (GIS) in each TMT batch, as described below. Similarly, GIS pooled standards were generated from the Banner, MSSB, Mayo, Aging, and UPenn cohorts. Samples were reduced with 1 mM dithiothreitol (DTT) at room temperature for 30 min, followed by 5 mM iodoacetamide (IAA) alkylation in the dark for another 30 min. Lysyl endopeptidase (Wako) at 1:100 (w/w) was added and digestion allowed to proceed overnight. Samples were then 7-fold diluted with 50 mM ammonium bicarbonate. Trypsin (Promega) was then added at 1:50 (w/w) and digestion was carried out for another 16 h. The peptide solutions were acidified to a final concentration of 1% (vol/vol) formic acid (FA) and 0.1% (vol/vol) trifluoroacetic acid (TFA), and desalted with a 30 mg HLB column (Oasis). Each HLB column was first rinsed with 1 mL of methanol, washed with 1 mL 50% (vol/vol) acetonitrile (ACN), and equilibrated with 2×1 mL 0.1% (vol/vol) TFA. The samples were then loaded onto the column and washed with 2×1 mL 0.1% (vol/vol) TFA. Elution was performed with 2 volumes of 0.5 mL 50% (vol/vol) ACN.

Mass Spectrometry Analysis for Label-free Proteomics

Mass spectrometry analyses of MSSB, ACT, BLSA, Banner, Mayo, and UPenn cohorts were performed on a Q-Exactive Plus mass spectrometer essentially as described¹⁷. Brain-derived tryptic peptides (2 μg) were resuspended in peptide loading buffer (0.1% FA, 0.03% TFA, 1% ACN) containing 0.2 pmol of isotopically labeled peptide calibrants (ThermoFisher 88321). Peptide mixtures were separated on a self-packed C18 (1.9 μm, Dr. Maisch, Germany) fused silica column (25 cm x 75 μM internal diameter; New Objective, Woburn, MA) by a NanoAcquity UHPLC (Waters, Milford, MA) and monitored on a Q-Exactive Plus mass spectrometer (ThermoFisher Scientific, San Jose, CA). Elution was performed over a 120 minute gradient at a rate of 400 nL/min with buffer B ranging from 3% to 80% (buffer A: 0.1% FA and 5% DMSO in water, buffer B: 0.1 % FA and 5% DMSO in ACN). The mass spectrometer cycle was programmed to collect one full MS scan followed by 10 data dependent MS/MS scans. The MS scans (300–1800 m/z range, 1,000,000 automatic gain control (AGC), 150 ms maximum ion time) were collected at a resolution of 70,000 at m/z 200 in profile mode, and the MS/MS spectra (2 m/z isolation width, 25% collision energy, 100,000 AGC target, 50 ms maximum ion time) were acquired at a resolution of 17,500 at m/z 200. Dynamic exclusion was set to exclude previous sequenced precursor ions for 30 seconds within a 10 ppm window. Precursor ions with +1 and +6 or higher charge states were excluded from sequencing.

Label-free Quantification

For the consensus LFQ search, 645 RAW files, including individual cases and pooled GIS samples from the MSSB, ACT, Banner and BLSA cohorts, were uploaded onto the Amazon Web Services (AWS) Cloud and analyzed using MaxQuant v1.6.3.4 with Thermo Foundation 2.0 for RAW file reading capability. The Mayo, BLSA precuneus, Aging, and UPenn cohorts were each searched separately using MaxQuant. The search engine Andromeda was used to build and search a concatenated target-decoy UniProt Knowledgebase (UniProtKB) containing both Swiss-Prot and TrEMBL human reference protein sequences (90,411 target sequences downloaded April 21, 2015), plus 245 contaminant proteins included as a parameter for the Andromeda search within MaxQuant. Methionine oxidation (+15.9949 Da), asparagine and glutamine deamidation (+0.9840 Da), and protein N-terminal acetylation (+42.0106 Da) were variable modifications (up to 5 allowed per peptide); cysteine was assigned a fixed carbamidomethyl modification (+57.0215 Da). Only fully tryptic peptides with up to 2 miscleavages were considered in the database search. A precursor mass tolerance of ±20 ppm was applied prior to mass accuracy calibration, and ±4.5 ppm after internal MaxQuant calibration. Other search settings included a maximum peptide mass of 6,000 Da, a minimum peptide length of 6 residues, and 0.05 Da tolerance for high resolution MS/MS scans. The false discovery rate (FDR) for peptide spectral matches, proteins, and site decoy fraction were each set to 1 percent. Quantification settings were as follows: re-quantify with a second peak-finding attempt after protein identification is complete; match full MS1 peaks between runs; use a 0.7 min retention time match window after an alignment function was found with a 20 minute retention time search space. The label free quantitation (LFQ) algorithm in MaxQuant^61,62 was used for protein quantitation. The quantitation method considered only razor and unique peptides for protein level quantitation. The total summed protein intensity was also used to assess overall signal drift across samples prior to LFQ normalization.

Isobaric Tandem Mass Tag (TMT) Peptide Labeling of ROS/MAP Brain Tissues

Prior to TMT labeling, cases were randomized by covariates (age, sex, PMI, diagnosis, etc.), into 50 total batches (8 cases per batch). Peptides from each individual case (n=400) and the GIS pooled standard (n=100) were labeled using the TMT 10-plex kit (ThermoFisher 90406). In each batch, TMT channels 126 and 131 were used to label GIS standards, while the 8 middle TMT channels were reserved for individual samples following randomization. Labeling was performed as previously described^21,37. Briefly, each sample (containing 100 μg of peptides) was re-suspended in 100 mM TEAB buffer (100 μL). The TMT labeling reagents were equilibrated to room temperature, and anhydrous ACN (256 μL) was added to each reagent channel. Each channel was gently vortexed for 5 min, and then 41 μL from each TMT channel was transferred to the peptide solutions and allowed to incubate for 1 h at room temperature. The reaction was quenched with 5% (vol/vol) hydroxylamine (8 μl) (Pierce). All 10 channels were then combined and dried by SpeedVac (LabConco) to approximately 150 μL and diluted with 1 mL of 0.1% (vol/vol) TFA, then acidified to a final concentration of 1% (vol/vol) FA and 0.1% (vol/vol) TFA. Peptides were desalted with a 200 mg C18 Sep-Pak column (Waters). Each Sep-Pak column was activated with 3 mL of methanol, washed with 3 mL of 50% (vol/vol) ACN, and equilibrated with 2×3 mL of 0.1% TFA. The samples were then loaded and each column was washed with 2×3 mL 0.1% (vol/vol) TFA, followed by 2 mL of 1% (vol/vol) FA. Elution was performed with 2 volumes of 1.5 mL 50% (vol/vol) ACN. The eluates were then dried to completeness using a SpeedVac.

High-pH Off-line Fractionation of ROS/MAP Brain Tissues

High pH fractionation was performed essentially as described⁶³ with slight modification. Dried samples were re-suspended in high pH loading buffer (0.07% vol/vol NH₄OH, 0.045% vol/vol FA, 2% vol/vol ACN) and loaded onto an Agilent ZORBAX 300 Extend-C18 column (2.1mm x 150 mm with 3.5 μm beads). An Agilent 1100 HPLC system was used to carry out the fractionation. Solvent A consisted of 0.0175% (vol/vol) NH₄OH, 0.01125% (vol/vol) FA, and 2% (vol/vol) ACN; solvent B consisted of 0.0175% (vol/vol) NH₄OH, 0.01125% (vol/vol) FA, and 90% (vol/vol) ACN. The sample elution was performed over a 58.6 min gradient with a flow rate of 0.4 mL/min. The gradient consisted of 100% solvent A for 2 min, then 0% to 12% solvent B over 6 min, then 12% to 40 % over 28 min, then 40% to 44% over 4 min, then 44% to 60% over 5 min, and then held constant at 60% solvent B for 13.6 min. A total of 96 individual equal volume fractions were collected across the gradient and subsequently pooled by concatenation⁶³ into 24 fractions and dried to completeness using a SpeedVac.

TMT Mass Spectrometry of ROS/MAP Brain Tissues

All fractions were resuspended in an equal volume of loading buffer (0.1% FA, 0.03% TFA, 1% ACN) and analyzed by liquid chromatography coupled to tandem mass spectrometry essentially as described⁶⁴, with slight modifications. Peptide eluents were separated on a self-packed C18 (1.9 μm, Dr. Maisch, Germany) fused silica column (25 cm × 75 μM internal diameter (ID); New Objective, Woburn, MA) by an Dionex UltiMate 3000 RSLCnano liquid chromatography system (ThermoFisher Scientific) and monitored on an Orbitrap Fusion mass spectrometer (ThermoFisher Scientific). Sample elution was performed over a 180 min gradient with flow rate at 225 nL/min. The gradient was from 3% to 7% buffer B over 5 min, then 7% to 30% over 140 min, then 30% to 60% over 5 min, then 60% to 99% over 2 min, then held constant at 99% solvent B for 8 min, and then back to 1% B for an additional 20 min to equilibrate the column. Buffer A was water with 0.1% (vol/vol) formic acid, and buffer B was 80% (vol/vol) acetonitrile in water with 0.1% (vol/vol) formic acid. The mass spectrometer was set to acquire in data dependent mode using the top speed workflow with a cycle time of 3 seconds. Each cycle consisted of 1 full scan followed by as many MS/MS (MS2) scans that could fit within the time window. The full scan (MS1) was performed with an m/z range of 350–1500 at 120,000 resolution (at 200 m/z) with AGC set at 4×10⁵ and maximum injection time 50 ms. The most intense ions were selected for higher energy collision-induced dissociation (HCD) at 38% collision energy with an isolation of 0.7 m/z, a resolution of 30,000, an AGC setting of 5×10⁴, and a maximum injection time of 100 ms. Five of the 50 TMT batches were run on the Orbitrap Fusion mass spectrometer using the SPS-MS3 method as previously described²¹.

TMT ROS/MAP Database Searches and Protein Quantification

All RAW files (1,200 RAW files generated from 50 TMT 10-plexes) were analyzed using the Proteome Discoverer suite (version 2.3, ThermoFisher Scientific). MS2 spectra were searched against the UniProtKB human proteome database containing both Swiss-Prot and TrEMBL human reference protein sequences (90,411 target sequences downloaded April 21, 2015), plus 245 contaminant proteins. The Sequest HT search engine was used and parameters were specified as follows: fully tryptic specificity, maximum of two missed cleavages, minimum peptide length of 6, fixed modifications for TMT tags on lysine residues and peptide N-termini (+229.162932 Da) and carbamidomethylation of cysteine residues (+57.02146 Da), variable modifications for oxidation of methionine residues (+15.99492 Da) and deamidation of asparagine and glutamine (+0.984 Da), precursor mass tolerance of 20 ppm, and a fragment mass tolerance of 0.05 Da for MS2 spectra collected in the Orbitrap (0.5 Da for the MS2 from the SPS-MS3 batches). Percolator was used to filter peptide spectral matches (PSMs) and peptides to a false discovery rate (FDR) of less than 1%. Following spectral assignment, peptides were assembled into proteins and were further filtered based on the combined probabilities of their constituent peptides to a final FDR of 1%. In cases of redundancy, shared peptides were assigned to the protein sequence in adherence with the principles of parsimony. Reporter ions were quantified from MS2 or MS3 scans using an integration tolerance of 20 ppm with the most confident centroid setting.

Parallel Reaction Monitoring (PRM) Analysis

Peptides from brain digests used for the first 3 batches of the untargeted UPenn cohort analysis (equal to 1 μg protein digestion) were used for targeted analysis on an Orbitrap Lumos™ Tribrid™ Mass Spectrometer (ThermoFisher Scientific) fitted with a Nanospray Flex ion source and coupled to a NanoAcuity liquid chromatography system (Waters). The tryptic peptides were resuspended in loading buffer (0.1% TFA, 500 ng/μl), and an external reference peptide mix (Promega) was spiked into the sample at the concentration of 0.5 pmol/μl. The solution (2 μl) was loaded onto a self-packed 1.9 μm ReproSil-Pur C18 (Dr. Maisch) analytical column (New Objective, 50 cm × 75 μm inner diameter; 360 μm outer diameter) heated to 60 °C. The capillary temperature and spray voltage was set at 300 °C and 2.0 kV, respectively. Elution was performed over a 100 min gradient at a rate of 350 nL/min with buffer B ranging from 1% to 32% (buffer A: 0.1% FA in water, buffer B: 0.1% FA in ACN). The column was then washed with 99% buffer B for 10 minutes and equilibrated with 1% B for 15 minutes. The mass spectrometer was set to collect in PRM mode using an inclusion peptide list (Supplementary Table 4). An additional full survey scan was collected to assess for possible interference. Full scans were collected at a resolution of 120,000 at 200 m/z with an AGC setting of 2×10⁵ ion and a maximum ion transfer (IT) time of 50 ms. For PRM scans, the settings were: resolution at 30,000 at 200 m/z, AGC target of 1×10⁵ ions, maximum IT time of 50 ms, microscans count of 1, isolation width of 1.6 m/z, and isolation offset of 0 m/z. A pre-optimized normalized collision energy of 32% was used to obtain the maximal recovery of target product ions. The top 5–10 product ions from this collision energy optimization were used for downstream peptide quantification.

Peptide Quantification

A spectral library was built using Skyline (version 4.2) based on tandem mass spectra gathered from previous data dependent acquisition methods. A Skyline template was then created to quantify the endogenous peptides. The template parameters were: centroided precursor mass analyzer, MS1 mass accuracy of 20 ppm; centroided product mass analyzer, MS/MS mass accuracy of 20 ppm; include all matching scans. All rawfiles were then imported and processed accordingly. The resulting extracted ion chromatograms (XICs) of selected fragments were manually inspected and peak picking adjustments were made accordingly. The sum of all product ion peak areas was calculated in Skyline and extracted for further statistical analyses. The peak areas were normalized using the peak areas of external reference peptides. Raw peptide intensities are provided at https://www.synapse.org/consensus.

Selected Reaction Monitoring (SRM) of ROS/MAP Brain Tissues

Samples were prepared for LC-SRM analysis using a standard protocol described elsewhere⁶⁵. Briefly, on average ~20 mg of DLPFC brain tissue from each subject was homogenized in denaturation buffer. After denaturation with DTT, 400 μg protein aliquots were taken for further alkylation with iodoacetamide followed by digestion with trypsin as described. The digests were cleaned using C₁₈ solid phase extraction, and 30 μL aliquots at 1 μg/μL concentrations were mixed with 30 μL synthetic peptide mix. LC-SRM experiments were performed on a nanoACQUITY UPLC (Waters) coupled to a TSQ Vantage mass spectrometer (ThermoScientific), with 2 μL of peptide injection for each brain sample. Buffer A was 0.1% FA in water and buffer B was 0.1% FA in 90% ACN. Peptide separations were performed on an Acquity UPLC BEH 1.7 μm C₁₈ column (75 μm i.d. × 25 cm) at a flow rate 350 nL/min using a gradient of 0.5% buffer B over 0 to 14.5 min, then 0.5% to 15% over 14.5 to 15.0 min, then 15% to 40% over 15 to 30 min, and then 45% to 90% B over 30 to 32 min. The heated capillary temperature and spray voltage was set at 350 °C and 2.4 kV, respectively. Both the Q1 and Q3 were set as 0.7 FWHM. A scan width of 0.002 m/z and a dwell time of 10 ms were used. All SRM data were analyzed using the Skyline software package. All data were manually inspected to ensure correct peak assignment and peak boundaries. The peak area ratios of endogenous light peptides and their heavy isotope-labeled internal standards (i.e., L/H peak area ratios) were then automatically calculated by the Skyline software, and the best transition without matrix interference was used for accurate quantification. Following homogenization of all tissues, small aliquots of protein from each of the samples was pooled, which were then digested and served as a global external pooled reference standard. Peptides generated from this pooled standard were scattered throughout the study (8 samples per 96-well plate) and were used to capture the technical variance that is due to sample preparation steps (except homogenization) and instrument measurements. The signal-to-noise ratio in quantification of each peptide was calculated as the ratio of variances across the human subject samples versus the technical controls. Peptides with a signal-to-noise ratio less than 2 were excluded from further analysis. The peptide relative abundances were log₂ transformed and centered at the median. The abundance of endogenous peptides was quantified as a ratio to spiked-in synthetic peptides containing stable heavy isotopes. The “light/heavy” ratios were log₂ transformed and shifted such that median log₂-ratio was zero. Normalization adjusted for differences in protein amounts among the samples. During normalization, the log₂-ratios were shifted for each sample to make sure the median was set at zero.

CSF Samples

All participants from whom CSF samples were collected provided informed consent under protocols approved by the Institutional Review Board (IRB) at Emory University. All patients received standardized cognitive assessments (including MoCA) in the Emory Cognitive Neurology clinic, the Emory Goizueta Alzheimer’s Disease Research Center (ADRC), and affiliated research studies (Emory Healthy Brain Study [EHBS] and Emory M²OVE-AD study). All diagnostic data were supplied by the ADRC and the Emory Cognitive Neurology Program. CSF was collected by lumbar puncture and banked according to 2014 ADC/NIA best practices guidelines. For patients recruited from the Emory Cognitive Neurology Clinic, CSF samples were sent to Athena Diagnostics and assayed for Aβ42, total-Tau, and phospho-Tau (CSF ADmark®) using the INNOTEST® assay platform. CSF samples collected from research participants in the ADRC, EHBS, and M²OVE-AD were assayed using the INNO-BIA AlzBio3 Luminex assay. In total, there were two cohorts of CSF samples that were used in the proteomics studies. Cohort 1 contained CSF samples from 150 healthy controls and 147 MCI/AD patients. Cohort 2 included CSF obtained from three groups: 32 cognitively normal, 31 AsymAD, and 33 MCI/AD. Cases and normal individuals with AsymAD were defined using established biomarker cutoff criteria for AD for each assay platform^66,67. Cohort information is provided in Supplementary Table 1.

CSF Protein Digestion

To generate peptides, all crude CSF samples were digested with LysC and trypsin. Briefly, 20 μL CSF from each sample was reduced and alkylated with 0.4 μL 0.5 M tris-2(-carboxyethyl)-phosphine (TCEP) and 2 μL 0.4 M chloroacetamide (CAA) with heating at 90°C for 10 min, followed by a 15 min water bath sonication. The samples were then further denatured by the addition of 67.2 μL of 8 M urea buffer (8 M urea, 100 mM NaHPO₄, pH 8.5) and digested overnight with 1.9 μg LysC (Wako) (1:10 enzyme to protein ratio according to the highest amount of sample). Following LysC digestion, the samples were diluted to 1 M urea using 50 mM ammonium bicarbonate. The same amount of trypsin (Promega) was then added (1:10 enzyme to protein ratio) and digestion was carried out for another 12 h. After trypsin digestion, the peptide solutions were acidified with a 1% TFA and 10% FA solution to a final concentration of 0.1% TFA and 1% FA. Peptides were desalted with a 30 mg C18 HLB column (Waters) and eluted in 1 mL of 50% ACN. Aliquots (120 μL) from cohort 1 (n=297) or cohort 2 (n=96) samples were pooled together and split into equal volume aliquots (880 μL) for use as the global internal standard (GIS) for TMT labeling. All samples and GIS were dried using a SpeedVac.

TMT Boost Channel

Signals of low abundant proteins in the TMT 11-plex were amplified using a boost channel, as previously described⁶⁸. A pooled CSF sample was created separately for each cohort by combining 50 μL from each sample in cohort 1 or cohort 2 into a pool for each cohort. Abundant proteins were removed using the High Select Top14 Abundant Proteins Depletion Resin (Thermo Scientific A36372BR) according to the manufacturer’s protocol, using a CSF-to-resin volume ratio of 1:1 and an incubation time of 15 min. After immunodepletion, protein concentrations were determined by BCA. Proteins were then reduced and alkylated (10 mM TCEP, 40 mM CAA) for 10 minutes at 90 °C. The samples were then subjected to bath sonication for 15 min and dried under vacuum in a SpeedVac. The immunodepleted pooled samples were re-suspended in 6 M urea buffer (6 M urea, 75 mM NaHPO₄, pH 8.5) at half the volume of the pooled sample prior to evaporation. Samples were digested overnight with LysC at an enzyme to protein ratio of 1:10. The following day, samples were diluted with 50 mM ammonium bicarbonate to reduce the urea concentration to 1M, and trypsin (Promega) was added (1:10 enzyme to protein ratio). Digestion was allowed to proceed for 12 hr. Peptides were then desalted using a 200 mg C18 Sep-Pak column, and the eluate was dried using a SpeedVac. Aliquots (600 μg) of the immunodepleted pooled CSF samples were separately dissolved in 100 mM TEAB buffer (625 μL) and labeled with 5 mg of TMT 126 channel reagent (cohort 1 lot# TF266326, cohort 2 lot# SG253447, ThermoFisher Scientific) in anhydrous ACN (256 μL). The reactions were allowed to proceed for 1 hr, and were subsequently quenched by adding 5% hydroxylamine (50 μL) and incubating for 15 min. The 126 channel was then added to the other channels, as described below.

TMT Labeling of Individual and GIS CSF Samples

All samples, including the GIS, were labeled with the 10-plex TMT kit plus an additional channel, for a total of 11 TMT channels (cohort 1 lot# TG273545 for 10-plex, TG273555 for channel 131C; cohort 2 lot# SI258088 for 10-plex, SJ258847 for channel 131C, ThermoFisher Scientific). Samples were grouped into batches as shown in Supplementary Table 1. The TMT labeling kit was equilibrated to room temperature and dissolved in anhydrous ACN (256 μL). The samples were reconstituted in 100 mM TEAB buffer (50 μL) and mixed with 0.4 mg (20.5 μL) of the corresponding labeling reagent. The labeling reactions were allowed to proceed for 1 hr, and were subsequently quenched with 5% hydroxylamine (4 μL). Per each TMT batch, labeled peptides from 9 channels (127N, 128N, 128C, 129N, 129C, 130N, 130C, 131, 131C) were mixed, desalted using a 100 mg C18 Sep-Pak column, and dried using a SpeedVac. The immunodepleted pooled sample labeled with the 126 channel (boost channel) was then added to each 9-channel TMT mixture at a ratio of 50:1 pooled to individual CSF sample by original volume:volume prior to evaporation. The sample mixtures were desalted using a 200 mg C18 Sep-Pak column, and dried using a SpeedVac.

Mass Spectrometry Analysis of CSF

All samples were resuspended in equal volume of loading buffer (0.1% FA, 0.03% TFA, 1% ACN). Peptide eluents were separated on a self-packed C18 (1.9 μm, Dr. Maisch, Germany) fused silica column (25 cm × 75 μM internal diameter (ID): New Objective, Woburn, MA) by an Easy-nLC system (ThermoFisher Scientific) and monitored on an Orbitrap Fusion Lumos mass spectrometer (ThermoFisher Scientific) interfaced with a high-field asymmetric waveform ion mobility spectrometry (FAIMS) Pro. Sample elution was performed over a 180 min gradient (buffer A: 0.1% FA in water, buffer B: 0.1% FA in 80% ACN) with flow rate at 225 nL/min. The gradient was from 1% to 8% buffer B over 3 min, then from 8% to 40% over 160 min, then from 40% to 99% over 10 min, and then held at 99% B for 10 min. The mass spectrometer was set to acquire data in positive ion mode using data dependent acquisition and three (−50, −65 and −85 V) different compensation voltages (CV). Data were acquired at each CV for 1 s during each cycle. Each cycle consisted of 1 full scan followed by as many MS2 and MS3 scans as possible within a 1 s timeframe. The full scan was performed with an m/z range of 450–1500 at 120,000 resolution (at 200 m/z) with an AGC setting of 4×10⁵ and maximum injection time 50 ms. The collision induced dissociation (CID) MS/MS scans were collected in the ion trap with an isolation window of 0.7 m/z, a collision energy of 35%, AGC setting of 1×10⁴, and a maximum injection time of 50 ms. The top 10 product ions were subjected to HCD synchronous precursor selection-based MS3 (SPS-MS3) as previously described²¹. For SPS-MS3 scans the isolation window was set to 2 m/z, the resolution to 50,000, the AGC to 1×10⁵, and the maximum injection time to 105 ms. For both cohorts, a single preliminary run of TMT batch 1 using the above parameters was used to create a target inclusion list of peptides that specifically excluded those from the top 15 most abundant proteins. This inclusion list was used for all TMT batches in cohort 1 (n=38) and in cohort 2 (n=12).

Database Searches and Protein Quantification of CSF

All RAW files were analyzed using the Proteome Discoverer Suite (version 2.3, ThermoFisher Scientific). MS/MS spectra were searched against the UniProtKB human proteome database (downloaded April 2015 with 90,411 total sequences). The Sequest HT search engine was used to search the RAW files, with search parameters specified as follows: fully tryptic specificity, maximum of two missed cleavages, minimum peptide length of 6, fixed modifications for TMT tags on lysine residues and peptide N-termini (+229.162932 Da) and carbamidomethylation of cysteine residues (+57.02146 Da), variable modifications for oxidation of methionine residues (+15.99492 Da), serine, threonine and tyrosine phosphorylation (+79.966 Da) and deamidation of asparagine and glutamine (+0.984 Da), precursor mass tolerance of 20 ppm, and a fragment mass tolerance of 0.6 Da. Percolator was used to filter PSMs and peptides to an FDR of less than 1%. Following spectral assignment, peptides were assembled into proteins and were further filtered based on the combined probabilities of their constituent peptides to a final FDR of 1%. In cases of redundancy, shared peptides were assigned to the protein sequence in adherence with the principles of parsimony. Reporter ions were quantified from MS3 scans using an integration tolerance of 20 ppm with the most confident centroid setting, as previously described²¹.

Controlling for Batch-specific Variance

We implemented a median polish algorithm for removing technical variance (e.g., due to tissue collection, cohort, or batch effects) from a two-way abundance-sample data table as originally described by Tukey⁶⁹. The algorithm is fully documented and available as an R function, which can be downloaded from https://github.com/edammer/TAMPOR. The algorithm implements iterations of the below equation, where batch and cohort are interchangeable.

Briefly, Equation 1 is applied to each protein measurement (LFQ or TMT reporter abundance) across all samples individually where the first term represents batch-wise median-centered abundance, and the second term is a batch-specific normalization factor comprised of the grand median of all batch-specific medians, divided by the appropriate batch-specific median of median-centered abundances. The data matrix is then log₂-transformed, and each log₂(ratio) is adjusted by subtraction of sample (column)-wise median log₂(ratio) for all proteins. Then, ratios are anti-logged and multiplied by the protein (row)-wise median of all samples used for the Eq. 1, term 1, denominator, extracted before Eq. 1 was executed. This process is iterated until convergence. The use of median polish ensures that the reduction of variance is robust to outliers while the overall algorithm preserves biological variance, given that batches have been randomized to avoid confounding batch with diagnosis or other biological traits. Prior to matrix assembly for the consensus analysis, intra-cohort batch effects were first removed in the MSSB (batch correction with 166 case samples across 7 batches) and Banner (batch correction with 178 case samples across 4 batches) cohorts. All remaining batch corrections restricted the first term denominator to global pooled (within cohort) standard sample abundances, and the second term used all individual case samples. Following removal of intra-cohort batch effects in MSSB and Banner, all samples were processed jointly with the algorithm in the same sample-protein matrix to capture biological variance across all samples in all four cohorts (ACT, Banner, BLSA, and MSSB) for the consensus analysis. The above algorithm was applied to a matrix in which proteins that had ≥ 50% missing values were removed. For the consensus LFQ network, 450 case samples (3 ACT outliers were removed prior to inclusion, as described below) classified as control, AsymAD, or AD by our unified criteria (see case classification methods above) were considered as “all samples” for denominators in Eq. 1. All remaining batch corrections listed as follows restricted the first term denominator to global pooled (within cohort) standard sample abundances, and the second term used all individual case samples.

For ROSMAP 50-batch TMT protein abundances, there were two pooled global internal standard channels in each TMT batch (n=64), and 400 individual case samples (non-internal standard samples). For the Hopkins aging cohort (84 case samples), global pool mixture samples (3 each per 3 batches) were used for the first term denominator, with the second term using all non-global pool mixture samples. For the UPenn PRM analysis (3 batches, 114 case samples, and 9 pooled controls), data were likewise batch corrected using 3 global pool mixture samples per batch for the first term denominator and all within-batch non-pooled samples for the second term. UPenn LFQ data (10 batches, 330 case samples, and 29 control pools) were similarly batch corrected as described. CSF 96-case and 300-case TMT normalized abundances were also batch corrected using the above algorithm, with equation 1 first term denominator restricted to global pooled (within cohort) standard sample abundances, while the second term used all individual non-internal standard case samples.

Regression for Covariates and Outlier Removal

No imputation of missing values was performed in any cohort. Nonparametric bootstrap regression was performed separately in each cohort by subtracting the trait of interest (age at death, sex, or postmortem interval (PMI)) times the median estimated coefficient from 1000 iterations of fitting for each protein in the cohort-specific log₂(abundance) matrix. Ages at death used for regression were uncensored. Case status/diagnosis was also explicitly modeled (i.e., protected) in each regression. Following regression of each individual cohort, we assessed whether any cohort-specific tissue dissection bias was present by performing a Spearman rank correlation of traits including age, sex, PMI, and white matter markers to the top five principle components (PC) of log₂(abundance). Network outlier case samples were not considered in the PCs, and were identified prior to PC analysis using Oldham’s ‘SampleNetworks’ v1.06 R script⁴ as previously published⁷⁰ using a 3 fold-SD cutoff of Z-transformed sample connectivity. The Spearman rank correlation was performed prior to correction of cohort-specific batch effects as described above, and after intra-cohort batch correction of the MSSB and Banner cohorts. All four of the cohorts were confirmed to have no significant PC correlation to age, sex, or PMI; however, ACT was observed to have a first PC significantly correlated (average rho=0.94) to protein abundance of white matter markers identified previously as oligodendrocyte coexpression network hubs¹⁷. These markers were BCAS1, SIRT2, MBP, and MAG. This white matter PC represented 27 percent of variance in the ACT cohort, whereas the white matter marker-correlated PC represented 7 to 12 percent variance in the other three cohorts. To adjust for this white matter variance in ACT, we applied a second round of bootstrap regression to the 62 non-outlier ACT case sample log₂(abundances), using the white matter PC as a regression covariate, and subtracted 28 percent of the white matter marker correlated variance to achieve a final variance of 12 percent after recalculation of the top 5 PCs. Abundance data for the 450 case samples were then assembled into a matrix of 3334 proteins, and cross-cohort batch correction by median polish was performed as described above. Finally, network outlier detection was performed as described above, which removed 31/450 cases from consideration in the four-cohort consensus network and differential abundance analyses. All outliers are listed in Supplementary Table 2. In all other cohorts that were not combined for the consensus network analysis, batch correction was performed first, followed by outlier removal, followed by removal of proteins with ≥ 50% missing values, and then regression of age, gender, and PMI prior to coexpression network and differential abundance analyses. In the Hopkins aging cohort, age was not considered as a trait for regression. In the CSF cohorts, only age at time of collection and sex were considered for regression.

Differential Expression Analysis

Differentially expressed proteins were found using one-way ANOVA followed by Tukey’s comparison post-hoc test across control, AsymAD and AD cases. Significantly altered proteins with corresponding p value are provided in Supplementary Tables 2 and 5 for consensus AD network proteins and astrocyte/microglial phenotype proteins, respectively. Differential expression is presented as volcano plots, which were generated with the ggplot2 package in R v3.5.2.

Protein Correlations to Aging

Protein expression levels for all 2,756 proteins measured in the aging cohort with fewer than 50% missing values were correlated to age at death using the bicor function, after regression for sex and PMI. In addition to bicor rho, the Student’s p value for significance of the correlation, FDR (q value), and signed z score for the correlation were calculated. The CSGene database of cellular senescence genes⁷¹ was cross-referenced with the protein data to annotate proteins involved in cellular senescence. This information is provided in Supplementary Table 3, along with differences between AD, control, or AsymAD for the same proteins in the consensus network, and their corresponding consensus network module relationships. Individual correlations for proteins in the consensus network to all traits provided in Supplementary Figure 1 are provided online at https://www.synapse.org/consensus.

Weighted Correlation Network Analysis (WGCNA)

We used the WeiGhted Correlation Network Analysis (WGCNA) algorithm for our network analysis pipeline, as previously described⁷⁰. A weighted protein co-expression network for the four-cohort consensus data was generated using the 3,334 log₂ protein abundance x 419 case-sample matrix that had undergone covariate and batch correction as described above. The WGCNA::blockwiseModules() function was used with the following settings: soft threshold power beta=5.5, deepSplit=4, minimum module size of 14, merge cut height of 0.07, mean TOM denominator, a signed network with partitioning about medioids (PAM) respecting the dendrogram, a minimum kME to remain in a module of 0.30, and a reassignment threshold of p<0.05. Specifically, we calculated pair-wise biweight mid-correlations (bicor, a robust correlation metric) between each protein pair, and transformed this correlation matrix into a signed adjacency matrix. The connection strength of components within this matrix was used to calculate a topological overlap matrix, which represents measurements of protein expression pattern similarity across cohort samples constructed on the pairwise correlations for all proteins within the network. Hierarchical protein correlation clustering analysis by this approach was conducted using 1-TOM, and initial module identifications were established using dynamic tree cutting as implemented in the WGCNA::blockwiseModules() function. Module eigenproteins were defined, which represent the most representative abundance value for a module and which explain covariance of all proteins within a module¹⁶. Pearson correlations between each protein and each module eigenprotein were performed; this module membership measure is defined as k_ME. After the initial network construction, 18 modules consisting of 14 or more proteins were detected. Given high k_ME similarity between some modules, we used the WGCNA::moduleMergeUsingKME() function to reduce the number of modules to 13, with the following parameters: percentage of module members checked for k_ME overlap of 50 percent (threshPercent=50), threshold for merging modules with high common k_{ME.intramodule} of 25 percent (mergePercent=25), and all other parameters default for the function. After module merging, MEs and signed kMEs were recalculated with the WGCNA::moduleEigengenes() and WGCNA::signedKME() functions, respectively. Finally, we ‘cleaned’ the network of assignments with aberrant k_ME (since WGCNA clustering into modules and merging of those modules use hybrid approaches not solely dependent on kME) by applying the following algorithm: remove module members with k_{ME.intramodule} <0.28, then check all grey (unclustered) proteins for any k_ME >0.35 and assign to the module with the highest k_{ME.intramodule}. For BLSA precuneus, Mayo temporal cortex, and ROS/MAP cohorts, we confirmed that they built networks using the pipeline described above, including batch correction, outlier removal, and age/sex/PMI regression. Network building was employed without module merging. The blockwiseModules function for ROS/MAP used power=10, minimum module size=30 (11,225 proteins in 323 non-excluded and non-outlier case samples); for BLSA precuneus used power=8, minimum module size=14 (3,348 proteins in 46 case samples); and for Mayo used power=5, minimum module size=14 (3,951 proteins in 107 case samples). For analysis of sex effects on the consensus network, the identical cases and analysis pipeline as described above were used to build the consensus network, except that covariance due to sex was not removed by regression prior to WGNCA analysis.

Network Preservation

We used the WGCNA::modulePreservation() function to assess network module preservation across cohorts. We also used this function to assess the effect of missing values on the consensus network. Z_summary composite preservation scores were obtained using the consensus network as the template versus each other cohort or missing value threshold tested, with 500 permutations. Random seed was set to 1 for reproducibility, and the quickCor option was set to 0. We also assessed network module preservation using synthetic eigenproteins. Briefly, protein module members in the consensus network template with a k_{ME.intramodule} among the top 20th percentile were assembled into a synthetic module in each target cohort, and synthetic modules with at least 4 members were used to calculate synthetic weighted eigengenes representing the variance of all members in the target network across case samples via the WGCNA::moduleEigengenes() function. Statistics and correlation scatterplots involving target cohort traits were then calculated and visualized.

Gene Ontology and Cell Type Marker Enrichment Analyses

To characterize differentially expressed proteins and co-expressed proteins based on gene ontology annotation, we used GO Elite v1.2.5 as previously published⁷⁰, with pruned output visualized using an in-house R script. Cell type enrichment was also investigated as previously published⁷⁰. Astrocyte subtype markers were obtained from Zamanian et al.²⁶. Microglia subtype markers were obtained from Rangaraju et al.²⁷.

GWAS Module Association

To determine if any protein products of GWAS targets were enriched in a particular module, we used the single nucleotide polymorphism (SNP) summary statistics from the International Genomics of Alzheimer’s Project (http://www.pasteur-lille.fr/en/recherche/u744/igap/igap_download.php) to calculate the gene level association value using MAGMA²⁵, as previously described¹⁷. APOE was added to the gene list and assigned a –log p value of 50, given its known strong association with AD. ApoE did not fall within a network module, and therefore did not influence the reported module enrichment results. Similar analyses were performed with GWAS candidates for Schizophrenia (SCZ) and Autism Spectrum Disorders (ASD). These GWAS datasets were provided and downloaded from the Psychiatric Genomics Consortium (http://www.med.unc.edu/pgc/downloads).

Curation of AD Mouse Model Purified Glial Transcriptomic and Proteomic Data

Orre et al.³² microarray data files (n=11; files last updated on July 20, 2015) were downloaded from the NCBI gene expression omnibus (GEO) datasets website (GEO identifier: GSE74615). The data, consisting of 22 samples, was normalized and centered based on limma R package user guidelines with backgroundCorrect method=”normexp” and normalizeBetweenArrays method=”quantile”. Transcripts with signal at least 110% of the 95^th percentile of all normalized red negative control spots on the arrays in at least 4/22 samples were considered as expressed above noise in at least one cell type (n=28,157 gene transcripts across N=14 purified microglial measurements and N=8 purified astrocyte measurements). WGCNA collapseRows function was used to handle multiple probe mappings to any gene, with probe selection method=”maxRowVariance”. Grubman et al.³³ peptide-level total peak area data for XO4⁺ (amyloid-β plaque phagocytosing) and XO4⁻ (amyloid-β plaque non-phagocytosing) acutely purified microglia from AD and WT mice, respectively, were summed and log₂-transformed to achieve protein-level data for 94 gene product proteins. Grubman et al. mRNA-level measurements were obtained directly and used without further processing.

Other Statistics

All statistical analyses were performed in R (v. 3.5.2). Boxplots represent the median, 25^th, and 75^th percentiles, and whiskers with staples represent measurements to the 5^th and 95^th percentiles of non-outlier samples; outlier samples are plotted beyond these whiskers as open circles. Correlations were performed using the biweight midcorrelation function as implemented in the WGCNA R package. Comparisons between two groups were performed by t test. Comparisons among three or more groups were performed with Kruskal-Wallis nonparametric ANOVA or standard ANOVA with Tukey or Dunnett post hoc pairwise comparison of significance. P values were adjusted for multiple comparisons by false discovery rate (FDR) correction where indicated. T-distributed Stochastic Neighbor Embedding (t-SNE) analysis was implemented with R package Rtsne, as previously published²⁷. Module membership graphs were generated using the igraph R package as previously described^17,72.