Large-scale proteome analysis of CSF implicates altered glucose metabolism in Alzheimer’s disease

Study participants

The data in this study came from two longitudinal AD cohorts of middle- and older-aged adults: the Wisconsin Registry for Alzheimer’s Prevention (WRAP)²⁶ and the Wisconsin Alzheimer’s Disease Research Center (ADRC)²⁷ (Table 1). Briefly, WRAP includes participants enriched for a parental history of AD dementia who were largely between the ages of 40 and 65 at the time of enrollment, fluent in English, able to perform neuropsychological testing, without a diagnosis or evidence of dementia at baseline, and without any health conditions that might prevent participation in the study. The ADRC study includes participants from one of several subgroups: mild late-onset AD, MCI, cognitively unimpaired middle-aged adults enriched for a parental history of presumed AD dementia, and age-matched healthy older controls (age > 65). Briefly, the ADRC participants were over the age of 45, with decisional capacity, and without a history of certain medical conditions (like congestive heart failure or major neurologic disorders other than dementia) or any contraindication to biomarker procedures. Participants in both the WRAP and ADRC cohorts were given diagnoses of AD, MCI, cognitively unimpaired, and others that were reviewed by a consensus review committee that included dementia-specialist physicians, neuropsychologists, and nurse practitioners²⁶. The National Institute of Neurological and Communicative Disorders and Stroke and Alzheimer’s Disease and Related Disorders Association (NINCDS-ADRDA)⁶ and NIA-AA⁷ criteria were used in defining the clinical diagnoses without reference to the participants’ CSF biomarker status. This study used the STROBE cohort reporting guidelines²⁸ and was performed as part of the GeneRations Of WRAP (GROW) study, which was approved by the University of Wisconsin Health Sciences Institutional Review Board. Participants in the ADRC and WRAP studies provided written informed consent.

CSF biomarkers

The CSF samples used for the biomarker analyses were acquired from lumbar punctures (LPs) using a uniform preanalytical protocol between 2010 and 2018 as previously described²⁹. Samples were collected in the morning using a Sprotte 24- or 25-gauge atraumatic spinal needle, and 22 mL of fluid was collected via gentle extraction into polypropylene syringes and combined into a single 30 mL polypropylene tube. After gentle mixing, samples were centrifuged to remove red blood cells and other debris. Then, 0.5 mL CSF was aliquoted into 1.5 mL polypropylene tubes and stored at -80 degrees Celsius within 30 minutes of collection.

All CSF samples were assayed between March 2019 and January 2020 at the Clinical Neurochemistry Laboratory at the University of Gothenburg. CSF biomarkers were assayed using the NeuroToolKit (NTK) (Roche Diagnostics International Ltd, Rotkreuz, Switzerland), a panel of automated Elecsys® and robust prototype immunoassays designed to generate reliable biomarker data that can be compared across cohorts. Measurements with the following immunoassays were performed on a cobas e 601 analyzer (Roche Diagnostics International Ltd, Rotkreuz, Switzerland): Elecsys β-amyloid (1–42) CSF (Aβ42), Elecsys Phospho-Tau (181P) CSF (ptau), and Elecsys Total-Tau CSF, β-amyloid (1–40) CSF (Aβ40), and interleukin-6 (IL-6). The remaining NTK panel was assayed on a cobas e 411 analyzer (Roche Diagnostics International Ltd, Rotkreuz, Switzerland), including markers of synaptic damage and neuronal degeneration (neurogranin, neurofilament light protein [NFL], and alpha-synuclein) and markers of glial activation (chitinase-3-like protein 1 [YKL-40] and soluble triggering receptor expressed on myeloid cells 2 [sTREM2]).

A total of nine established CSF biomarkers for AD were analyzed in this study: the Aβ42/Aβ40 ratio, ptau, the ptau/Aβ42 ratio, NFL, alpha-synuclein, neurogranin, YKL-40, sTREM2, and IL-6. Since the CSF biomarker measurements were to be used as outcomes, each biomarker was assessed for skewness using the skewness function of the R package moments (version 0.14)³⁰. Any biomarker with a skewness ≥ 2 was transformed with a log₁₀-transformation to better meet the normality assumption of regression. The outcomes that were log₁₀-transformed were NFL and IL-6.

Samples used in this study were then assigned to pathological categories from the NIA-AA ATN research framework¹¹ using binary cut-offs for CSF amyloid and tau positivity. The development of these research cut-offs is described in detail elsewhere²⁹. Briefly, cut-offs were estimated via ROC analysis on a subsample of n = 185 participants (cognitively impaired and unimpaired) who underwent [11C] PiB-PET imaging within two years of an LP. Using the Matlab perfcurve function³¹ with an equally weighted cost function³², the optimal Aβ42/Aβ40 threshold was 0.046 and the optimal ptau/Aβ42 threshold was 0.038. Thresholds for ptau181 were determined by establishing a reference group of 223 CSF amyloid (Aβ42/Aβ40) negative, cognitively unimpaired younger participants (ages 40-60 years). Biomarker positivity thresholds for these analytes were set at +2 standard deviations (SD) above the mean of this reference group (ptau threshold = 24.8 pg/mL). In this study, A+ and T+ were defined based on the CSF Aβ42/Aβ40 and ptau thresholds, respectively. The final pathological categories for this study included amyloid negative and tau negative (A-T-); amyloid positive and tau negative (A+T-); and amyloid positive tau positive (A+T+). The fourth possible category of amyloid negative and tau positive (A-T+) was not included in this study as these samples are considered to represent non-AD pathological change¹¹.

CSF metabolomics

All samples used in this study had CSF metabolomics data available from the WRAP and ADRC cohorts. The details of the CSF sample collection, handling, and metabolomics profiling have been previously described^{33, 34}. Briefly, fasting CSF samples were drawn from study participants in the morning through LP and then mixed, centrifuged, aliquoted, and stored at -80 degrees Celsius. Samples were kept frozen until they were shipped overnight to Metabolon, Inc. (Durham, NC), which similarly kept samples frozen at -80 degrees Celsius until analysis. Metabolon used Ultrahigh Performance Liquid Chromatography-Tandem Mass Spectrometry (UPLC-MS/MS) to conduct an untargeted metabolomics analysis of the CSF samples. The metabolites were then annotated with metabolite identifiers, chemical properties, and pathway information. Metabolite measurements were divided by the median measurement for that metabolite across all samples. Missing values for xenobiotic metabolites were imputed to 0.0001, while missing values for non-xenobiotic metabolites were imputed to half of the minimum value among all other measured samples for that metabolite.

The initial data set contained 412 metabolites from 1,172 CSF samples across 687 unique individuals. A total of 13 metabolites that were missing for ≥ 50% of samples were removed. One sample was removed for missing ≥ 40% of metabolite values. A total of 9 metabolites with low variance (interquartile range = 0) were then removed. A log₁₀ transformation was applied to all metabolite values. A total of 220 samples from a clinical trial were excluded from analysis. The processed data set contained 390 metabolites quantified on 951 CSF samples from 609 unique individuals, including all of the individuals here on whom proteomics data were generated.

Genome-wide genotyping

Genome-wide genotypes were also available for all samples in this study. The genotyping in both the WRAP and ADRC had been previously conducted³⁵. For the WRAP cohort, DNA from whole blood samples were genotyped with the Illumina Multi-Ethnic Genotyping Array at the University of Wisconsin Biotechnology Center³⁴. Pre-imputation quality control (QC) steps included removing samples and variants with a high missingness (> 5%) or inconsistent genetic and self-reported sex. Samples from individuals of European descent were imputed using the Michigan Imputation Server³⁶ and the Haplotype Reference Consortium (HRC) reference panel³⁷. Variants with poor quality (R² < 0.8) or out of Hardy-Weinberg equilibrium (HWE) were removed after imputation, leaving a total of 1,198 samples with 10,499,994 single nucleotide polymorphisms (SNPs). In the ADRC, whole blood samples were genotyped by the Alzheimer’s Disease Genetics Consortium (ADGC) at the National Alzheimer’s Coordinating Center (NACC) using the Illumina HumanOmniExpress-12v1_A, Infinium HumanOmniExpressExome-8 v1-2a, or Infinium Global Screening Array v1-0 (GSAMD-24v1-0_20011747_A1) BeadChip assay. Initial quality control was conducted on each chip’s data separately, removing variants or samples with high missingness (> 2%), out of HWE (P < 1x10^-6), or with inconsistent genetic and self-reported sex. The remaining samples were then imputed with the Michigan Imputation Server, phased using Eagle2³⁸, and imputed to the HRC reference panel. As before, variants of low quality (R² < 0.8) or out of HWE were removed. The data sets from the different chips were then merged together, leaving a data set with 377 samples of European descent and 7,049,703 SNPs. The WRAP and ADRC data sets were then harmonized to each other and to the 1000 Genomes Utah residents with Northern and Western European ancestry (CEU)³⁹ data set, using the GRCh37 genome build. Ambiguous SNPs were removed, and then the remaining SNPs were aligned to the same strand and allele orientations as the ADRC data set.

The 137 samples from this study were then extracted from this combined genetic data set and further processed using PLINK⁴⁰ (v1.90b6.3). To ensure sufficient data were available for use in the prediction models, only SNPs with no missing data and with a minor allele count of 20 or greater among the 137 samples were retained. Linkage disequilibrium (LD) pruning was then applied using a window size of 1000 kb, an R² threshold of 0.1, and the 1000 Genomes CEU samples as the reference data set. The pruning resulted in a data set of 38,652 SNPs.

APOE genotyping

Each sample was additionally assigned an APOE genotype based on the participant’s combination of the ε2, ε3, and ε4 alleles for APOE from a separate set of genotyping. DNA was extracted from whole blood samples, which was then genotyped for the APOE alleles using competitive allele-specific PCR-based KASP genotyping for rs429358 and rs7412³³.

Proteomics sample selection

Based on the results of our pilot study for CSF proteomics²⁵, we had estimated a priori that a sample of approximately 150 would be sufficient for 80% power to detect most of the observed protein-AD diagnosis associations from the pilot using the R package pwr (version 1.3-0)⁴¹ (Supplementary Figure 1). The process of selecting samples for CSF proteomics generation began by considering all CSF samples from fasted, successful LPs (n = 1,440) from 823 unique participants across WRAP and ADRC. From there, each CSF sample was matched to its closest set of CSF biomarker data, CSF metabolomics data, and consensus conference diagnosis. Samples were excluded if there was insufficient material for proteomics analysis, if they were part of a clinical trial, or if they had been used already in our pilot study. To simplify the downstream analyses, only one sample (the most recent) per participant was considered when there were multiple samples. An approximately equal number of samples per AT-defined subgroup (A-T-, A+T-, A+T+) was selected, prioritizing samples with available genomic data and metabolomic data. A total of 140 samples were selected to have proteomics data generated.

Protein extraction and digestion

CSF protein concentration was determined by protein BCA assay (Thermo Scientific). CSF aliquots were moved to 96-plates and dried down using a SpeedVac Concentrator (Thermo Scientific) before being resuspended in a lysis buffer consisting of 10 mM TCEP, 40 mM CAA, 100 mM Tris pH 8, and 8M urea. The sample solution was then diluted to 25% strength using 100 mM Tris pH 8 before the addition of protease. Trypsin was added to the protein solution at an approximate ratio of 50:1 w/w and digested overnight at ambient temperature. The digestion reaction was quenched by acidification using TFA. Digested peptides were desalted using Strata-X Polymeric Reverse Phase plates (Phenomenex) before being dried down in the SpeedVac Concentrator overnight. Dried down samples were resuspended in 0.2% FA and peptide concentration was determined using a peptide BCA (Thermo Scientific). Peptide samples were injected directly from the 96-well plates.

Offline fractionation

Pooled samples for each of the three disease groups were created by combining 10 µL of CSF from each sample in that disease group. These three pooled samples were then prepared using the extraction and digestion protocol described above. The three desalted, digested peptide solutions were then fractionated using high-pH reverse-phase liquid chromatography. Separation was performed using an Agilent Infinity 2000 HPLC with a 150 mm C18 reverse-phase column (Waters, XBridge Peptide BEH, particle size 3.5 µm). Mobile phase buffer A was a freshly prepared mixture 10 mM ammonium formate pH 9.5, and mobile phase buffer B was a freshly prepared mixture of 80% MeOH, 10 mM ammonium formate pH 9.5. The gradient method was 20 minutes in length with fractions collected from minute 5 to minute 20, with a flow a rate of 800 nL/min across the entire method. The method initiated with a concentration of 5% B before increasing to 35% by minute 2. Percent B increased to 100% by 13 minutes. From 5 to 20 minutes, 32 fractions were collected in round-bottom 96 well plates in a time-based manner. Fractions were concatenated into a total of 16 by combining every other column in the collection plate. Fractionated samples were injected directly from the collection plate.

Online chromatography

We used a method previously developed in our pilot study²⁵: a single-shot nano-liquid chromatography-tandem mass spectrometry (nLC-MS/MS) method for quantitative and fast analysis of CSF protein extracts. Reverse phase columns used online with the mass spectrometer were packed using an in-house column packing apparatus described previously⁴². In brief, 1.5 µm BEH particles were packed into a fused silica capillary purchased from New Objective (PicoTip, Stock # PF360-75-10-N-5) at 30,000 psi. During online LC separations, capillary was heated to 50° C and interfaced with mass spectrometer via an embedded emitter. For online chromatography, a Dionex UltiMate 300 nanoflow UHPLC was used with mobile phase A consisting of 0.2% FA and mobile phase B consisted of 70% ACN, 0.2% FA. A flow rate of 310 nL/min was used throughout with the method increasing from 0% to 7% B over the first four minutes. Percent B then increased to 49% B by 59 minutes before a wash step of 100% B from 62 to 67 minutes. Method finished with an equilibration step from minute 68 to 78 of 0% B.

Tandem mass spectrometry

Peptides eluting from the column were ionized by electrospray ionization and analyzed using a Thermo Orbitrap Eclipse hybrid mass spectrometer. Survey scans were collected in the Orbitrap at a resolution of 240,000 with a normalized AGC target of 250% (1e6) with Advanced Precursor Determination engaged across the range of 300–1400 m/z. Precursors were isolated for tandem MS scans using a window of 0.5 m/z, with a dynamic exclusion duration of 22 seconds and a mass tolerance of 15 ppm. Precursors were dissociated using HCD with a normalized collision energy of 25%. Tandem scans were taken over the range 130–1350 m/z using the “rapid” setting with a normalized AGC target of 300% (3e4) and a maximum injection time of 18 ms.

The resulting raw data files were searched in MaxQuant^{43, 44} using fast LFQ and a full human proteome with isoforms downloaded from UniProt (downloaded June 14, 2017). Oxidation of methionine and acetylation of the N terminus were allowed as variable modifications, and carbamidomethylation of cysteine was set as a fixed modification. Proteins were searched using an FDR of 1% with a minimum peptide length of 7 and a 0.5 Da MS/MS match tolerance. Matching between runs was utilized, applied with a retention time window of 0.7 minutes. Protein abundance data were extracted in the form of LFQ Intensity from the “proteinGroups.txt” output file. Throughout this manuscript, each protein group is referred to by the first listed majority protein from its annotation from MaxQuant. The protein data were annotated with Entrez IDs (via R package org.Hs.eg.db⁴⁵, version 3.11.4), UniProt⁴⁶ IDs, and gene information (GENCODE⁴⁷, version 37, and the HUGO Gene Nomenclature Committee, HGNC, database⁴⁸). When the gene annotations conflicted or were absent from one of these databases for a given UniProt ID, the gene identifiers were taken in the order of resources listed.

Proteomics quality control

After removing several samples with injection or other technical issues, the proteomics data set included 2,040 proteins across 137 samples. These data underwent a strict quality control pipeline: proteins that were missing for 33% or more of samples (either overall or within an AT group) were removed; samples missing 33% or more of proteins were removed; and proteins with an interquartile range of 0 were removed (Supplementary Figure 2, Supplementary Figure 3). A total of 137 samples with 915 proteins remained (Supplementary Table 1). The label-free quantification (LFQ) values for each protein were then log₂-transformed. The remaining missing values were then randomly imputed based on a normal distribution derived from the lower end of the observed values for that protein (the observed distribution mean was shifted by -1.8 and the SD shrunk by a factor of 0.3) (Supplementary Figure 4). This imputation was performed separately within each AT group. Finally, each protein was Z-score transformed.

Proteomics descriptive analysis

The set of proteins quantified in the CSF was extensively profiled (Supplementary Table 2). The pairwise correlation of all proteins was calculated (nominally significant results with correlation P < 0.05 shown in Supplementary Table 3) and then visualized with a heatmap with hierarchical clustering to show the underlying patterns of covariation (R package ComplexHeatmap⁴⁹, version 2.4.3) (Figure 1a). The structure was further examined with a principal components analysis (PCA), scree plot (Supplementary Figure 5), and plot of the first two PCs by AT group (R package factoextra⁵⁰, version 1.0.7) (Supplementary Figure 6) to assess the presence of independent signals among the CSF proteome and their relationship to the AT groups. A pathway analysis was then performed to examine the differences between the set of proteins quantified in the CSF and the overall human proteome. The enrichment of Gene Ontology (GO) terms⁵¹, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways^{51, 52}, and Disease Ontology (DO) gene sets⁵³ among the CSF proteome against the entire human proteome was calculated using the R packages clusterProfiler⁵⁴ (version 3.16.1), DOSE⁵⁵ (version 3.14.0), and ReactomePA⁵⁶ (version 1.32.0) (Figure 1b, Figure 1c, Supplementary Table 4). The presence of clusters was then assessed with Gaussian mixture modeling using the R package mclust⁵⁷ (version 5.4.6). The number of clusters (3) was chosen based on the elbow of the plot of the Bayesian Information Criterion (BIC) (Supplementary Figure 7). Enrichment of gene set ontologies across the clusters was repeated with the GO, KEGG, and DO sets and plotted (Figure 1d, Supplementary Table 5).

• Download figure
• Open in new tab

Figure 1: CSF proteomics descriptive analyses

a) A correlation heatmap of the CSF proteins (n = 915) across all participants (n = 137) is shown. The dendrogram above the heatmap shows the results of hierarchical clustering of the proteins. An intricate set of correlation patterns can be seen, with both large clusters of proteins (e.g., top-left and bottom-right of the plot) and small local clusters seen throughout. b) An enrichment plot of the GO biological process terms among quantified proteins in the CSF compared to the entire human genome is shown. Significantly enriched processes included extracellular processes, processes involving axons and synapses, and immune system processes. c) A network representation of enriched DO terms among the CSF proteome is shown. The tan nodes represent significantly enriched disease terms, and the gray nodes represent proteins whose genes were associated with those terms. Three clusters of terms emerged: Alzheimer’s disease and tauopathy; coronary artery disease; and arteriosclerotic cardiovascular disease and arteriosclerosis. d) A cluster comparison plot of KEGG pathways between the three main clusters of CSF proteins is shown. The clusters were generated by Gaussian mixture modelling. Cluster one shows an enrichment of extracellular matrix and cell cycle pathways; cluster two shows enrichment of immune system-related pathways; and cluster three shows enrichment of a handful of other pathways.

Protein-AT category associations

The primary objective was the identification of differentially expressed proteins across the three AT groups. This analysis was performed using an analysis of covariance (ANCOVA) model comparing each protein across the three groups, controlling for age at LP and sex (Supplementary Table 6). A Bonferroni correction for the number of proteins tested (P = 0.05 / 915 = 5.46 × 10^-5) was used for reporting significant results. The distributions of the top-associated proteins across the AT spectrum were plotted (Figure 2a). To assess whether signal enrichment was likely due to an artifact, the ANCOVA analyses were repeated with randomly permuted AT group labels. A quantile-quantile (Q-Q) plot was generated to assess the presence of signal enrichment across the proteome for AT-related differences and to compare the permuted and non-permuted analyses (Figure 2b). Since the APOE gene is known to have a significant effect on AD risk, we examined whether APOE genotype was driving the observed AT-protein associations. The ANCOVA analyses were repeated but with the count of APOE ε4 alleles included as an additional covariate. The same Bonferroni correction was used as before. The resulting AT-associated proteins were compared to the results from the original ANCOVA analyses (Supplementary Figure 8). The set of associated proteins from the non-permuted analysis was then assessed for enriched GO, KEGG, and DO gene sets against the human proteome as before (Figure 2c, Supplementary Table 7).

• Download figure
• Open in new tab

Figure 2: Proteins associated with AT group

a) The distributions of the top 10 most significantly associated proteins with AT group are shown. A number of different patterns were seen across increasing AT pathology, including proteins that increased consistently, decreased consistently, or increased only with tau positivity. b) The quantile-quantile (Q-Q) plots of the protein-AT ANCOVA association tests are shown above, with the distribution of P values shown separately for the original (“Not permuted”) and permuted data sets. Substantial signal enrichment was seen across the CSF proteome, with that signal absent in the permuted data set. c) The enriched KEGG pathways among the AT-associated proteins are shown, revealing a general perturbation of metabolic pathways.

To examine the direction of effect of each protein, a logistic regression was performed with A+T+ (vs. A-T-) as the outcome and a protein as the main predictor, controlling for age at LP and sex and using the same Bonferroni threshold for significance as the ANCOVA analyses. The sample size for the logistic regression was smaller (n = 98) due to the exclusion of the A+T-samples. The overlap between the set of significantly associated proteins and the set of significantly associated proteins from the ANCOVA analysis was displayed in a Venn diagram (R package ggVennDiagram⁵⁸, version 1.0.7) (Supplementary Figure 9), and the odds ratios were presented in a volcano plot (Supplementary Figure 10).

Protein-CSF biomarker associations

Each protein was then tested for association with each of the 9 CSF biomarkers (the Aβ42/Aβ40 ratio, ptau, the ptau/Aβ42 ratio, NFL, alpha-synuclein, neurogranin, YKL-40, sTREM2, and IL-6). Linear regression models were used to regress each CSF biomarker on each protein, controlling for age at CSF sample and sex and using a Bonferroni correction for the total number of tests (9 × 915 = 8,234 tests; P = 0.05 / 8,234 = 6.07 × 10^-6) (Table 2, Supplementary Table 8). The results were summarized with a Q-Q plot showing the signal enrichment for each biomarker along with a sensitivity analysis where the regression models were repeated with the biomarker values randomly permuted to test the robustness of each biomarker’s signal enrichment (Supplementary Figure 11). The cross-biomarker relationships among the significantly associated proteins were then visualized as a bipartite graph using the R package tidygraph (version 1.2.0)⁵⁹ and the Fruchterman-Reingold algorithm. A community structure network analysis was performed using the greedy hierarchical agglomeration algorithm⁶⁰ implemented in igraph (version 1.2.5)⁶¹ to identify clusters of proteins among the protein-biomarker associations. An upset plot was then created showing the set of significantly associated proteins unique to each subset of biomarkers (Supplementary Figure 12) using the UpSetR package (version 1.4.0)⁶². Pathway enrichment analyses were performed as before comparing each biomarker’s set of significantly associated proteins (Supplementary Table 9).

Replication of top results in the Knight ADRC

A replication data set from the Knight ADRC was used to validate findings from the main analyses. The Knight ADRC data set included samples from the brain (n = 360 participants), CSF (n = 717), and plasma (n = 490). The recruited individuals from the Knight ADRC cohort were evaluated by Clinical Core personnel of Washington University. For brain samples, a brain autopsy was performed by the Knight ADRC neuropathological core, and AD status was determined by a postmortem neuropathological analysis. Neuropathological phenotypes including Braak tau, Consortium to Establish a Registry for Alzheimer’s Disease (CERAD) Aβ, α-synuclein pathology, postmortem interval (PMI), age at onset (AAO), age at death, and brain weight were obtained for all brain samples. For individuals with CSF and plasma data, cases received a clinical diagnosis of AD in accordance with standard criteria, and AD severity was determined using the clinical dementia rating (CDR) scale⁶³ at the time of lumbar puncture (for CSF samples) or blood draw (for plasma samples). Brain tissues were collected from fresh frozen human parietal lobes. CSF and plasma samples were collected in the morning after an overnight fast and immediately centrifuged and stored at -80°C until assayed⁶⁴. The Institutional Review Boards of Washington University School of Medicine in St. Louis approved the study; research was performed in accordance with the approved protocols.

For deep omics characterization in the brain, CSF, and plasma tissues, the level of 1,305 proteins were quantified using the SOMAscan assay, a multiplexed, aptamer-based platform⁶⁵. The assay covers a dynamic range of 10⁸ and measures all three major categories: secreted, membrane, and intracellular proteins. The proteins cover a wide range of molecular functions and include proteins known to be relevant to human disease. As previously described by Gold et al.⁶⁵, modified single-stranded DNA aptamers are used to bind specific protein targets that are then quantified by a DNA microarray. Protein concentrations are quantified as relative fluorescent units (RFU). Aliquots of 150 μL of tissue were sent to the Genome Technology Access Center at Washington University in St. Louis for protein measurement.

Quality control was performed at the sample and aptamer levels using control aptamers (positive and negative controls) and calibrator samples. At the sample level, hybridization controls on each plate were used to correct for systematic variability in hybridization. The median signal over all aptamers was used to correct for within-run technical variability. This median signal was assigned to different dilution sets within each tissue. For brain and CSF samples, a 20% dilution rate was used. For plasma samples, three different dilution sets (40%, 1%, and 0.005%) were used.

As described in detail⁶⁶, additional quality control was performed by identifying and removing protein and analyte outliers by applying four criteria: 1) Minimum detection filtering. If the analyte for a given sample was less than the limit of detection (LOD), the sample was deemed an outlier. Collectively, if the number of outliers given an analyte was less than 15% of the total sample size, the analyte was kept. 2) Flagging analytes based on the scale factor difference. 3) Coefficient of variation (CV) of calibrators lower than 0.15, where the CV for each aptamer was calculated as the standard deviation divided by the mean of each calibrator at the raw protein level. 4) Interquartile range (IQR) strategy. Outliers were identified if the subject was located 1.5-fold of the IQR outside of either end of the distribution given the log₁₀-transformation of the protein level. Analytes were kept after passing all the criteria above for the downstream statistical analysis. An orthogonal approach was used to call subject outliers based on IQR. After this second removal of analytes, subject outliers were examined and removed again.

To obtain the proteomic signatures of sporadic AD status, differential abundance analysis was performed by using linear regression of the log-transformed protein levels. In all three tissues, sporadic AD status was considered as a main predictor. In each tissue, we performed surrogate variable analysis while protecting status and age to correct for unmeasured heterogeneity⁶⁷. Age at death or at measurement, gender, and the resulting surrogate variables were included as covariates. To validate identified proteins using the related traits, analyses using AAO and AD neuropathology characteristics (Braak scores and CDR at death) were considered for brain data, AAO and CSF ptau/Aβ ratio were considered for CSF data, and AAO was considered for plasma data.

The minimum principal component (PC) number that cumulatively explained 95% of the variance for each tissue after QC was obtained. The number of PCs was 169, 230, and 75 in CSF, plasma, and brain data, respectively. These numbers corresponded to 2.96 × 10^-4, 2.17 × 10^-4, and 6.67 × 10^-4, respectively, as the Bonferroni-corrected P value thresholds. The use of these thresholds was more conservative than the false discovery rate (FDR).

The Knight ADRC analyses were used as a replication data set for the top findings from the main analyses performed in the Wisconsin ADRC and WRAP, focusing on the significantly associated proteins from the top implicated biological pathway from the protein-AT group and protein-biomarker analyses. The association of these proteins were compared to the results from the Knight ADRC association tests conducted in the CSF, brain, and plasma to see if their associations and directions of effect were replicated using the significance thresholds defined above (Table 3, Supplementary Table 10).

Secondary analysis of insulin-related proteins

Based on the results of the AT category and biomarker associations and the pathway analysis suggesting a relationship with glucose regulation (described below), the set of proteins excluded during the QC process due to low sample size was examined for proteins related to insulin signaling pathways, including any of the GLUT proteins (SLC2A family), insulin (INS), insulin receptor (INSR), insulin-like growth factor 1 (IGF1), IGF-1 receptor (IGF1R), insulin receptor substrate 1 (IRS1), IRS 2 (IRS2), phosphoinositide 3-kinase (PI3K), RAC-alpha serine/threonine-protein kinase (AKT1), mechanistic target of rapamycin (mTOR), and glycogen synthase kinase 3 (GSK3A). Proteins that failed the missingness threshold of 33% but were present for 50% or more of samples were investigated further but without the use of imputed data points. The relationships between the proteins and AT category (Supplementary Figure 13) and the CSF biomarkers (Supplementary Figure 14) were plotted, with ANCOVA and linear regression analyses to test for association between the proteins and AT group and the CSF biomarkers performed as previously (Supplementary Table 11).

• Download figure
• Open in new tab

Figure 3: Network analysis of protein-biomarker associations

A bipartite graph representation of the proteins significantly associated with the CSF biomarkers after Bonferroni correction is shown. The nodes representing the biomarkers are larger and in green. Proteins are represented as smaller nodes, with an edge to a biomarker representing a significant association between a protein and a biomarker. The colors of the protein nodes and underlying shaded regions correspond to three distinct communities identified with the fast greedy modularity optimization algorithm. Three such communities were identified. Community 1 included proteins associated with the more traditional AD biomarkers, such as ptau, ptau/Aβ42, and Aβ42/Aβ40. Community 2 centered around neurogranin, including many proteins uniquely associated with neurogranin. Community 3 included proteins associated with the remaining biomarkers of alpha-synuclein, YKL-40, NFL, and sTREM2. No proteins were significantly associated with IL-6. Many proteins were shared in the center of the network, particularly among neurogranin, ptau, and alpha-synuclein.

• Download figure
• Open in new tab

Figure 4: Multiomic predictive model performance for amyloid and tau

The ability of different omic data sets to predict the key AD CSF biomarkers for amyloid and tau is summarized. a) The receiver operator characteristic (ROC) curves for Aβ42/Aβ40, ptau/Aβ42, and ptau are shown for the held-out testing data set predictions. Each omic prediction set is plotted with a different line. b-d) A 2D density plot summarizes the actual values and predicted probability of positivity for each CSF biomarker (Aβ42/Aβ40, ptau/Aβ42, and ptau, respectively) from the proteomic predictor models applied to the held-out testing data set. The text labels in the corners refer to the prediction categories: TN = true negative, FP = false positive, TP = true positive, and FN = false negative. The vertical red line indicates the threshold for a hard classification of biomarker positivity from the prediction models. The horizontal red line indicates the binary threshold for the CSF biomarker determined from previous work (see Methods).

Multiomic prediction of amyloid and tau

The CSF proteomic data set was then combined with the CSF metabolomic, genome-wide genotyping, and demographic (age at sample and sex) data sets. After the quality control steps described previously for each ome, all 137 samples had values for all of the multiomic features (915 proteins, 390 metabolites, 38,652 SNPs, and 2 demographic features). The multiomic data set was then used to predict different biomarker positivity states²⁹: Aβ42/Aβ40-positive, ptau-positive, and ptau/Aβ42-positive. Each ome (CSF proteome, CSF metabolome, genome, and demographics) was used individually along with a fifth multiomics predictor set (comprising all omes) to predict each outcome with an elastic net⁶⁸ model (R package glmnet⁶⁹, version 4.0-2; alpha parameter = 0.5).

For each biomarker and predictor pair, the procedure was the same. First, one-third of the data was held out as a testing set and the remaining two-thirds used as the training set. Within the training set data, 100-iteration, 3-fold cross-validation was used to select the best lambda value (11 possible values ranging from 10^-5 to 1) according to AUC using the tidymodels⁷⁰ R package (version 0.1.3). The best-performing model was then run on the entire training data set using the chosen lambda and used to predict the outcome on the held-out testing data set. The performances of the different omic models were then compared with ROC curves and 2D histograms showing the raw biomarker levels against the predicted classifications for each biomarker for each subject (Figure 4). The mean model metrics across each of the 4,000 folds were calculated (Supplementary Table 12).

Sample summary

CSF samples from 137 WRAP and ADRC participants were selected as described in the Methods, roughly evenly distributed across the three AT groups of interest (Table 1, Supplementary Figure 1). Most (102, 74.5%) of the participants were cognitively unimpaired at the time of the sample, with 16 (11.7%) and 19 (13.9%) participants having an MCI or AD dementia diagnosis, respectively. The age and sex distributions across the AT categories varied, with worse AT pathology having a higher average participant age and a greater proportion of males. The amyloid and tau measures reflected the AT categorizations as expected. The remaining CSF biomarkers showed a general increase with increasing AT pathology with the exception of IL-6, which fluctuated across the groups.

CSF proteomics descriptive analyses

The nLC-MS/MS analysis, MaxQuant identification, and LFQ quantification generated a total of 2,040 protein groups across the participants. After the proteomics quality control steps (Supplementary Table 1, Supplementary Figure 2, Supplementary Figure 3, Supplementary Figure 4), 915 proteins remained (Supplementary Table 2). Included in these proteins were YKL-40 (correlation with immunoassay measurement = 0.352, P = 2.40 × 10^-5), sTREM2 (correlation with immunoassay measurement of sTREM2 = 0.490, P = 1.26 × 10^-9), apolipoprotein E (APOE), amyloid precursor protein (APP; correlation with Aβ42 = 0.136, P = 0.114), amyloid-like protein-1 (APLP1), and APLP2. The tau protein was not reliably quantified by nLC-MS/MS in our samples. Little difference in protein missingness was seen by AT group (Supplementary Figure 2, Supplementary Figure 3). The CSF proteome showed a rich correlation structure with both larger clusters and smaller pockets of highly correlated proteins (Figure 1a, Supplementary Table 3). Further interrogation with PCA underscored this complexity, with the first 4 PCs collectively explaining only half (49.89%) of the total variance (Supplementary Figure 5), with the top 2 PCs not explained by either AT or sex (Supplementary Figure 6).

Pathway enrichment analysis comparing the proteins quantified in the CSF to the entire human proteome revealed significant enrichment of terms related to extracellular, neuronal, immune system, and platelet pathways (Figure 1b, Supplementary Table 4). Significantly enriched DO pathways among the cohort included three groups of proteins with some small overlap: Alzheimer’s disease and tauopathy; coronary artery disease; and arteriosclerotic cardiovascular disease and arteriosclerosis (Figure 1c).

Given the apparent presence of clusters of proteins based on the correlation structure, the CSF proteome was divided into 3 clusters based on a Gaussian mixture model (Supplementary Figure 7). These three clusters were then compared to each other for the differential enrichment of biological pathways (Supplementary Table 5). The KEGG terms revealed a pattern where the smallest cluster (2) was enriched for immune system and cholesterol pathways while the other two larger clusters were enriched for extracellular and metabolism-related pathways (Figure 1d).

Protein-AT associations

The ANCOVA tests revealed 61 statistically significant associations between proteins and AT group after multiple testing correction (P < 5.46 × 10^-5), with a total of 496 (54.2%) of the proteins nominally associated (P < 0.05) (Supplementary Table 6). The differences in distribution of the top ten proteins revealed a number of different patterns in relation to amyloid and tau pathology (Figure 2a). Some proteins increased (S4R371, SMOC1) or decreased (FBLN1) consistently as AT pathology increased. Other proteins did not change from A- to A+ but did change from T- to T+ (AATC, ALDOA, GUAD, DDAH1, CRYM, and F5H5Q2). Overall, there was an enrichment of statistical signal across the proteome that was not seen in the permutation sensitivity analysis (Figure 2b). Controlling for the APOE ε4 allele count did not substantially change the results, with 53 of the 61 proteins remaining significantly associated when the APOE variable was added to the ANCOVA models (Supplementary Figure 8). Among the set of significantly associated proteins, a number of biological pathways were enriched relative to the full human proteome (Supplementary Table 7), including extracellular matrix, secretory granule, and vesicle lumen GO cellular component terms; peptidase regulation GO molecular function terms; and glucose metabolic KEGG pathways (Figure 2c).

When the logistic regression model was used to test for the direction of effect between proteins and A+T+ (vs. A-T-), only 9 proteins were significantly associated with being A+T+, and all of these proteins were also significantly associated in the ANCOVA model (Supplementary Figure 9). All but 1 (FBLN1) of the proteins significantly associated with A+T+ were increased in A+T+ relative to A-T-(Supplementary Figure 10).

Protein-CSF biomarker associations

When each of the 9 CSF biomarkers (Aβ42/Aβ40, ptau, ptau/Aβ42, NFL, alpha-synuclein, neurogranin, YKL-40, sTREM2, and IL-6) was regressed on each CSF protein, a total of 636 protein-biomarker associations were statistically significant after Bonferroni correction (P < 6.07 × 10^-6; Supplementary Table 8). As with the protein-AT associations, there was widespread association signal across the proteome with the CSF biomarkers that was not seen in the permutation test, except for IL-6, which had no significantly associated proteins (Supplementary Figure 11). The top 3 significantly associated proteins per biomarker are summarized in Table 2. A total of 119 significantly enriched pathways among the biomarker-specific sets of significantly associated proteins was observed, with glucose metabolic pathways noted to be enriched among amyloid-related biomarkers (Supplementary Table 9).

The network plot and subsequent community analysis of the protein-biomarker associations revealed three communities (modularity = 0.256) among the network (Figure 3). One community largely comprised the more traditional AD biomarkers of ptau, ptau/Aβ42, and Aβ42/Aβ40; a second community centered around the proteins associated with neurogranin; and the third community included the remaining biomarkers of alpha-synuclein, YKL-40, NFL, and sTREM2 (IL-6 had no significant protein associations and was not included). The largest number of shared associations across the biomarkers occurred among neurogranin, ptau, and alpha-synuclein, which shared 103 protein associations between at least two of those biomarkers (Supplementary Figure 12).

Replication of top pathway results in the Knight ADRC

The glucose metabolism pathway (REACTOME ID R-HSA-70326) was significantly enriched among all three of the amyloid and tau measures, so the results of the 9 proteins from this pathway that were significantly associated with one of the amyloid or tau biomarkers were chosen for replication in the Knight ADRC. These proteins included MDH1, ALDOA, PGK1, TPI1, PGAM1, PKM, GOT1, ALDOC, and ENO1. Of these 9 proteins, 4 of them (MDH1, ALDOA, PGK1, TPI1) were replicated in the Knight ADRC CSF case-control analysis, with MDH1 (P = 1.07 × 10^-19) and ALDOA (P = 7.43 × 10^-14) both meeting the Knight ADRC’s Bonferroni-corrected threshold for significance (Table 3). The directions of effect were concordant except for PGK1. There was no significant association of these proteins with AD case status in the brain or plasma samples. Among the other AD-related outcomes in the Knight ADRC data, MDH1 and ALDOA were both statistically significantly associated with the CSF ptau/Aβ ratio (P = 8.08 × 10^-28 and P = 1.76 × 10^-27, respectively), matching what was seen for these proteins with this biomarker in the Wisconsin cohorts both in terms of strength of association (P = 3.21 × 10^-10 and P = 7.34 × 10^-8, respectively) and direction of effect (positive) (Supplementary Table 10).

Secondary analysis of insulin-related proteins

Of the insulin-related proteins of interest, only IGF-1 and AKT1 were identified in the proteomics workflow. AKT1 was only quantified in one subject and was thus not suitable for further analysis, but IGF-1 was quantified in 82 samples (59.9%) and analyzed further using only the non-imputed measurements. A trend in missing values by AT category was noted: 44.6% of samples were missing IGF-1 in A-T-, 41.0% in A+T-, and 33.3% in A+T+. The ANCOVA analysis of IGF-1 did not show a statistically significant difference of the protein across AT categories (P = 0.170), though the distribution of the protein appeared to increase with amyloid positivity (Supplementary Figure 13). The association analysis between IGF-1 and the CSF biomarkers revealed a nominally significant negative association with Aβ42/Aβ40 (P = 0.011) and positive association with ptau/Aβ42 (P = 0.009) (Supplementary Table 11, Supplementary Figure 14).

Multiomic prediction models for amyloid and tau

The results of the multiomic amyloid and tau prediction models revealed a consistent pattern where the CSF proteome outperformed the other omic data sets in predicting positivity based on the core biomarkers of Aβ42/Aβ40, ptau/Aβ42, and ptau (Figure 4, Supplementary Table 12). The predictive model based on the CSF proteome (number of predictors selected ranged from 37-97) achieved a high AUC across for all three biomarkers (Aβ42/Aβ40 AUC = 0.924, ptau/Aβ42 AUC = 0.917, and ptau AUC = 0.891), performing slightly better than even the integrative model. For Aβ42/Aβ40 and ptau positivity, the sensitivity and specificity values further demonstrated the relative superiority of the proteomics model. For Aβ42/Aβ40 positivity, the sensitivity of the CSF proteome model was 0.947 compared to much lower values (0.316–0.579) from the other models. Similarly, for ptau positivity, the specificity of the CSF proteome model (0.571) was much higher than the other models (0.000–0.143). In all cases, the genome-based model performed poorly (AUCs ranged from 0.500–0.628) (Figure 4a). The 2D histograms showing the performance of the CSF proteome relative to the raw biomarker values highlighted the effective classification by the proteomic models with effective delineation between positive and negative amyloid statuses (Figure 4b–d).