The genetics of circulating BDNF: towards understanding the role of BDNF in brain structure and function in middle and old ages.

Brain-derived neurotrophic factor (BDNF) plays an important role in brain development and function. Substantial amounts of BDNF are present in peripheral blood, and may serve as biomarkers for Alzheimer's disease incidence as well as targets for intervention to reduce Alzheimer's disease risk. With the exception of the genetic polymorphism in the BDNF gene, Val66Met, which has been extensively studied with regard to neurodegenerative diseases, the genetic variation that influences circulating BDNF levels is unknown. We aimed to explore the genetic determinants of circulating BDNF levels in order to clarify its mechanistic involvement in brain structure and function and Alzheimer's disease pathophysiology in middle-aged and old adults. Thus, we conducted a meta-analysis of genome-wide association study of circulating BDNF in 11 785 middle- and old-aged individuals of European ancestry from the Age, Gene/Environment Susceptibility-Reykjavik Study (AGES), the Framingham Heart Study (FHS), the Rotterdam Study and the Study of Health in Pomerania (SHIP-Trend). Furthermore, we performed functional annotation analysis and related the genetic polymorphism influencing circulating BDNF to common Alzheimer's disease pathologies from brain autopsies. Mendelian randomization was conducted to examine the possible causal role of circulating BDNF levels with various phenotypes including cognitive function, stroke, diabetes, cardiovascular disease, physical activity and diet patterns. Gene interaction networks analysis was also performed. The estimated heritability of BDNF levels was 30% (standard error = 0.0246, P-value = 4 × 10-48). We identified seven novel independent loci mapped near the BDNF gene and in BRD3, CSRNP1, KDELC2, RUNX1 (two single-nucleotide polymorphisms) and BDNF-AS. The expression of BDNF was associated with neurofibrillary tangles in brain tissues from the Religious Orders Study and Rush Memory and Aging Project (ROSMAP). Seven additional genes (ACAT1, ATM, NPAT, WDR48, TTC21A, SCN114 and COX7B) were identified through expression and protein quantitative trait loci analyses. Mendelian randomization analyses indicated a potential causal role of BDNF in cardioembolism. Lastly, Ingenuity Pathway Analysis placed circulating BDNF levels in four major networks. Our study provides novel insights into genes and molecular pathways associated with circulating BDNF levels and highlights the possible involvement of plaque instability as an underlying mechanism linking BDNF with brain neurodegeneration. These findings provide a foundation for a better understanding of BDNF regulation and function in the context of brain aging and neurodegenerative pathophysiology.

Introduction

Brain-derived neurotrophic factor (BDNF), a member of the neurotrophins superfamily, has received remarkable attention in Alzheimer’s disease research and is thought to play a key role in its pathophysiology. Extensive research in animal models suggests that BDNF plays central roles in normal brain functions and development, regulates growth and survival of neurons, promotes long-term potentiation and modulates synaptic transmission and activity-dependent plasticity (Kokaia ; Murer ; Lee ; Chan ; Waterhouse ). BDNF and its receptor, tyrosine receptor kinase B, are highly expressed and widely distributed throughout the CNS, especially in the cerebral cortex and hippocampus (Phillips ; Wetmore ) where they are implicated in memory formation and long-term potentiation (Alonso ; Bekinschtein ). Furthermore, brain BDNF expression has been reported to modify the relation of Alzheimer’s disease pathology to cognitive decline (Buchman ).

Substantial amounts of BDNF are also present in peripheral blood and are highly correlated with CSF levels, as 75% of BDNF is produced in neurons and glia (Rasmussen ), and readily crosses the blood–brain barrier (Pan ; Karege ). Thus, there is a growing interest in circulating BDNF as a marker of Alzheimer’s disease, which may be of high applicability due to easy access to blood samples. Indeed, evidence demonstrates reduced levels of circulating BDNF in persons with Alzheimer’s disease and mild cognitive impairment (Song ), as well as in healthy individuals who are destined to develop Alzheimer’s disease (i.e. during Alzheimer’s disease prodromal phase) (Weinstein ). Furthermore, decreased BDNF levels in human serum have been associated with poorer cognitive function, smaller hippocampal volume and widespread brain amyloid burden (Hwang ).

Besides its value as a potential Alzheimer’s disease biomarker, BDNF is thought to mediate the association between lifestyle and Alzheimer’s disease risk, because it is induced by factors such as reduced caloric intake (Lee ) and increased physical activity (Szuhany ). This potential role is of high importance in light of recent efforts to ameliorate Alzheimer’s disease risk by conducting healthier lifestyle, which is estimated to account for ∼30% of Alzheimer’s disease cases worldwide (Baumgart ). Lastly, BDNF may potentially have therapeutic effects in Alzheimer’s disease, which can be reached through lifestyle changes or by administration of exogenous BDNF or stimulation of its receptors (Tapia-Arancibia ; Lee ).

Despite the plethora of evidence suggesting a fundamental role of circulating BDNF in Alzheimer’s disease, little is known about its genetic determinants. An exception is the functional polymorphism (rs6265) in the BDNF gene that results in a valine-to-methionine substitution at codon 66 (Val66Met) which has been extensively studied with regard to neurodegenerative disease (Shen ), however, it is not clear if it affects the expression of BDNF in peripheral blood (Terracciano ; Kaess ). Hence, we carried out a genome-wide association study (GWAS) and explored heritability and post-GWAS functional annotations in the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) cohorts to better understand the genetic determinants of circulating BDNF levels.

Materials and Methods

Study population

We performed a GWAS meta-analysis of 11 785 participants of European ancestry from four cohorts: Age, Gene/Environment Susceptibility-Reykjavik Study (AGES), Framingham Heart Study (FHS), Rotterdam Study and Study of Health in Pomerania (SHIP-Trend). All are part of the CHARGE consortium. A detailed description of each study can be found in the Supplementary material. Written informed consent was obtained from all participants. Each study was approved by local ethical committees or the institutional review boards (Supplemental material for details).

Definition of phenotypes

BDNF levels were measured in serum in all studies but the Rotterdam Study in which levels were measured from plasma. Information on BDNF source and assay type is provided in Supplementary Table 1.

Genotyping and imputation

Genotyping was performed using a variety of commercial arrays across the participating studies. Information on genotyping platforms, quality control procedures and imputations methods for each participating study are provided in Supplementary Table 2. Similar quality control procedures were applied for each study (Supplementary Table 2). Each study performed genotype imputations using Haplotype Reference Consortium (HRC, version 1.1) or 1000 Genomes (phase 1, version 3) imputed genotype dosages. A total of 9 950 208 single-nucleotide polymorphisms (SNPs) were included in the study.

Genome-wide association study

Each participating study performed the GWAS of circulating BDNF levels under an additive model using variant allele dosage as predictors and circulating BDNF concentrations as dependent variable. The association analyses were adjusted for age, sex, familial relationship (in family-based studies) or study site (in multi-site studies). Population stratification was controlled for by including principal components derived from genome-wide genotype data. Study-specific details on covariates are provided in Supplementary Table 3. The effect sizes of each cohort were then combined using the fixed-effect effective-sample-size weighted Z-score meta-analysis in METAL, which take into account possible differences in BDNF sources and assays (Willer ). Effective-sample-size was defined as the multiple of sample size and observed divided by expected variance for imputed allele dosage. Variants with minor allele frequency (MAF) <0.01 were filtered out. After the meta-analysis, 250 variants that surpassed the P-value < 5×10−8 threshold were included in the clumping analysis by PLINK (http://zzz.bwh.harvard.edu/plink, 20 October 2020, date last accessed) (Purcell ) to identify the independent top SNPs. The clumping procedure started with the most significant SNP as the index SNP, formed a clump of all SNPs that had r2 ≥ 0.1 with index SNP and also within 500 kb distance from the index SNP. This process was then iteratively applied to the rest of the SNPs not already clumped until all SNPs were clumped. The index SNPs as well as all other SNPs within its 500 kb region were illustrated in the regional association and recombination plots, using LocusZoom (Pruim ) based on the hg 19 University of California, Santa Cruze (UCSC) Genome Browser assembly for visualization of candidate genes and variants in linkage disequilibrium (LD). Based on the GWAS summary data, we then performed multi-SNP-based conditional and joint association analysis using Genome-wide Complex Trait Analysis - Conditional & Joint (GCTA-COJO) (Yang ). This procedure performs a stepwise model selection to select independently associated SNPs. P-value and collinearity cut-offs were set to 5×10−8 and 0.9, respectively, and the distance for complete linkage equilibrium was set to 10 Mb. Because individual-level genotype data were not available, we used the FHS genotype data as the reference sample.

We also conducted a fixed-effect meta-analysis on SNP effect sizes with unit being one standard deviation of BDNF levels in each cohort using inverse-variance weighted method in METAL to obtain an overall effect estimate, which was then compared with effect sizes from individual cohort. The previously implied rs6265 (VAL66MET) was among the SNPs in the GWAS. Lastly, in a sensitivity analysis, we restricted our GWAS to serum samples only (excluding the Rotterdam Study in which BDNF levels were measured from plasma) since plasma and serum BDNF levels may indicate distinct biological relevance (Gejl ).

Polygenic heritability

Heritability of circulating BDNF levels was estimated in the FHS, a population-based cohort with family structure. We used sequential oligogenic linkage analysis routines to estimate the total heritability, i.e. inter-person variability in the phenotype explained by the proportion of shared genetics between the persons. Specifically, a linear mixed-effects model was fitted with age and sex as fixed effects and with person-specific random effects. Between families, these random effects are uncorrelated. Within a family, the correlation of the random effects between each pair of subjects is twice their kinship coefficient. The total heritability is the variance of the random effects divided by the total phenotypic variance. S equential oligogenic linkage analysis routines do not use any genotype data in the heritability estimation but require individual-level data and knowledge of family relationships. The assumption is of no residual correlations caused by other environmental or genetic factors. There are other methods [e.g. genomic restricted maximum likelihood and linkage disequilibrium score regressions (Yang ; Bulik-Sullivan )], which estimate the so-called narrow sense heritability in which the heritability is the variance explained by all SNPs in the GWAS. When most variants are common as in this study, narrow sense heritability is expected to be lower than the total heritability because rare causal variants are typically poorly tagged by common variants.

Local SNP heritability

The single SNP heritability estimation was estimated using beta estimation and allele frequency based on formula: . The proportion of variance explained by all top SNPs together was also estimated using Heritability Estimation from Summary Statistics (HESS) (version 0.5.3- beta, https://huwenboshi.github.io/hess, 20 October 2020, date last accessed) based on GWAS summary data while accounting for LD among variants (Petyuk ). Fixed-effect effective-sample-size weighted meta-analysed Z-scores were taken as GWAS summary association input. Genome partition for each chromosome was conducted to create approximately LD-independent loci. The genome partition file was provided by Berisa and Pickrell (2016). We additionally included the non-overlapping 60 kb region of the identified top independent SNPs in the genome partition file in order to calculate local SNP heritability of top hits. 1000 Genomes for Europeans was used as reference panel. All SNPs in the reference panel have MAF greater than 1%.

Functional annotations

Expression quantitative trait loci in peripheral blood mononuclear cell, monocyte, brain and T cells

To determine whether the genome-wide significant variants influence gene expression, we searched if those variants were also expression quantitative trait loci (eQTL) using data of peripheral blood mononuclear (PBMC), monocyte, brain and T cells from multiple sources and projects. Functionally annotated SNPs were mapped to genes based on physical position on the genome and eQTL associations (all tissues). To correct for multiple testing, the false discovery rate (FDR) <0.05 was used as a threshold.

T cells and monocytes

In 211 subjects from the Brigham and Women’s Hospital PhenoGenetic Project, peripheral venous blood was obtained from healthy control volunteers. Flow-sorted CD14+CD16− monocyte expression data collected using the Affymetrix GeneChip Human Gene 1.0 ST Array and genotype data imputed to MAF >0.01 from healthy subjects of EA. The raw expression intensity values were normalized using RMA normalization. More detailed methods on data collection have been described previously (Shi ).

Peripheral blood mononuclear cell

Three hundred and sixty-three subjects with demyelinating disease and frozen PBMC samples were selected from the CLIMB (Raj ) at the Partners Multiple Sclerosis Center in Boston. Of these subjects, 232 were paired with imputed genotype data, and eQTL analysis was performed. Analysis pipeline of PBMCs mirrored that of T cells and monocyte eQTL analysis as previously described. Collection of expression data was previously described (Gauthier ).

Dorsolateral prefrontal cortex

In 407 subjects from the two longitudinal cohorts of aging persons, the Religious Orders Study and the Rush Memory and Aging Project (Bennett ; De Jager ). eQTL analysis was performed on bulk tissue from the dorsolateral prefrontal cortex. RNA was extracted from the grey matter of dorsolateral prefrontal cortex, and next-generation RNA sequencing (RNA-Seq) was done on the Illumina HiSeq for samples with an RNA integrity score >5 and a quantity threshold >5 μg. Expression was quantile-normalized by fragments per kilobase of transcript per million fragments mapped (FPKM), correcting for batch effect with Combat. These adjusted FPKM values were used for analysis. Genotyping was performed on either the Affymetrix GeneChip 6.0 platform (1878 participants, 909 600 SNPs) or the Illumina OmniQuad Express platform (456 participants, 730 525 SNPs). DNA was extracted from whole blood, lymphocytes, or frozen brain tissue, as previously described. Details of RNA-seq and genotype collection were previously described (De Jager ).

Module eQTL analysis

In a set of 508 subjects from the Religious Orders Study (ROS) and the Rush Memory and Aging Project (MAP) (Bennett ; De Jager ), the SpeakEasy (SE) clustering algorithm was used to define 47 gene expression modules, or groups of co-expressed genes (Mostafavi ). Normalized gene expression data were generated from bulk tissue dissected at autopsy from the dorsolateral prefrontal cortex. The modules range in size from between 20 and 556 gene members (median = 331). The mean normalized expression level of the genes within each module was calculated to generate a single meta-feature per module. In 494 subjects with both genetic and expression data, GWAS was then performed on each of the 47 module meta-features. The generation of the modules and RNA-seq has been previously described (Gaiteri ). To correct for multiple testing, an FDR <0.1 was used as threshold.

Associations of SNP with neuropathological phenotypes (neuritic plaque, neurofibrillary tangles, amyloid and tau) in ROS and MAP

We related the following phenotypes to the lead SNPs influencing circulating BDNF levels in brains from ROS and MAP samples: (i) neuritic plaques, (ii) neurofibrillary tangles, (iii) β-amyloid load and (iv) tau tangle density (n = 1100). To correct for multiple testing, an FDR<0.1 was used as a threshold.

ROS and MAP are prospective cohort studies of aging where study participants who are cognitively normal at enrolment agree to annual clinical evaluations and to Anatomic Gift Act donating their brains at the time of death. They provide an informed consent as well as a repository consent which allows the data to be repurposed. Each subject undergoes a detailed quantitative neuropathological examination, detailed ante-mortem clinical and neuropsychological profiling, and banking of PBMCs. The follow-up rate of survivors exceeds 90% and the autopsy rate exceeds 80%. A further description can be found elsewhere (Bennett ; De Jager ).

Brain autopsies were performed and each brain was inspected for common pathologies relating to loss of cognition in aging populations (Bennett ; Ottoboni ). Neurofibrillary tangles and neuritic plaques were visualized by modified Bielschowsky silver stain, then counted and scaled in five brain regions. Composite scores for each of these pathology types were derived by scaling the counts within each of the five regions, and taking the square root of the average of the regional scaled values to account for their positively skewed distribution (Bennett ; Ottoboni ). β-amyloid load and tau tangle density were measured in eight brain regions by immunohistochemistry, a composite score was calculated, which was square root transformed (Replogle ).

Protein quantitative trait loci in ROS and MAP proteomics

In 920 subjects from the ROS and MAP (Bennett ; De Jager ), Selection Reaction Monitoring (SRM) proteomics was performed using frozen tissue from dorsolateral prefrontal cortex (Yu ). In all, 228 peptides were quantified and analysed, corresponding to 148 proteins, selected via a curated list of genes implicated in neurodegenerative disease. The samples were prepared for LC-SRM analysis using standard protocol. All the SRM data were analysed by Skyline software. All the data were manually inspected to ensure correct peak assignment and peak boundaries. The peak area ratios of endogenous light peptides and their heavy isotope-labelled 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. The peptide relative abundances were log2 transformed and centred at the median. GWAS was performed on all of the 228 peptide measurements. To correct for multiple testing, an FDR<0.1 was applied.

Associations of gene expression with neuropathological phenotypes with in ROS and MAP

We further related neuropathological findings (definition detailed in part c) to the expression of the 14 gene set (gene within 60 kb of GWAS lead variants influencing circulating BDNF levels) in 508 brains from the Religious Orders Study and Rush Memory and Aging Project (ROSMAP) samples. With 5 phenotypes and 14 genes, a Bonferroni corrected P-value of 0.05/(5×14) = 0.000714 was applied to control for multiple testing, and maintain an alpha level of 0.05.

Gene prioritization and pathway analysis

To establish functional connections, we conducted three different analyses implemented in the DEPICT (v1, https://data.broadinstitute.org/mpg/depict, 20 October 2020, date last accessed) (Pers ). First, to prioritize genes with relevant biological roles in the circulating BDNF-associated loci, we tested functional similarities among genes from different associated regions where genes with high functional similarity across regions obtained lower prioritization P-values. Second, we analysed expression enrichment across particular tissues or cell types by testing whether genes in the associated circulating BDNF loci had high expression in any of the 209 Medical Subject Headings (MeSH) annotations, using data from 37 427 expression arrays. Third, we performed a gene set enrichment analysis to test whether the genes in the associated circulating BDNF loci were enriched in reconstituted gene sets. The 10 968 gene sets tested were generated from diverse databases, including Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome KnowledgeBase (REACTOME), the InWeb database (high-confidence protein–protein interaction), and the Mouse Genetics Initiative (phenotype–genotype relationships). In all three analyses, we used the FDR to adjust for multiple testing; significance was defined at FDR=0.20. The DEPICT analyses were based on independent lead SNPs (r2 < 0.3) with P-values < 1×10−5.

Gene interaction network analysis

We took 22 unique genes located within 60 kb of the genome-wide significant SNPs as input to Ingenuity Pathway Analysis system to generate networks using information in Ingenuity Knowledge Base curated from the published literature consisting of thousands of genes and gene products that interact with each other. Interactions between all pairs of genes or gene products in the Ingenuity Knowledge Base have been abstracted into the Global Molecular Network with genes or gene products as nodes and the interactions between them as edges connecting the nodes. Using Ingenuity Pathway Analysis system, the 22 focus genes of BDNF were first sorted based on their triangle connectivity (the number of triangles in the Global Molecular Network that contains the gene). Starting from the top-ranked focus gene, neighbourhood focus genes were added to form a network, optimizing the specific connectivity (how much a new added focus gene’s neighbourhood overlaps the current network) until the maximum network size of 35 was reached or no other focus gene was in the neighbourhood of the network. The top-ranked focus gene was removed from the focus genes set and the previous process was repeated for the remaining focus genes to form networks. After networks were constructed for each focus gene, from the smallest constructed network, a small connected focus gene network (number of node ≤17) was combined to another small network through a linker gene (the linker gene was selected so that it had maximal edges to both networks), until the merged network’s size was larger than 17 or no linker gene exited to connect the network to other small networks. Additional genes were then added to the periphery of the network to provide additional biological context to those focus genes, based on the number of the edges that connect them to any gene in the constructed network. Finally, the networks were ranked by P-scores, calculated as -log10 (P-value) of the Fisher’s exact test. In a constructed network with n genes and f focus genes, the P-value was the probability of getting f or more focus genes in a set of n genes randomly selected from the Global Molecular Network.

Mendelian randomization

We performed Mendelian randomization (MR) analyses using the genome-wide loci identified as instruments to test the hypothesis that circulation BDNF levels are causally associated with cardiovascular diseases, stroke, general cognitive function, Alzheimer’s disease, smoking, depression, diabetes, blood lipid levels and brain MRI. In addition, we conducted MR analyses to test whether circulating BDNF levels are causally related to physical activity (Doherty ) and diet patterns (Cole ). Of note, despite the strong correlation between global cognition and intelligence phenotypes, participants’ consent did not allow the study of intelligence in FHS and other CHARGE cohorts. Therefore, we did not include genetic studies of intelligence in our MR evaluation of circulating BDNF levels and global cognition. Two-sample MR analyses were performed using the R package MR (https://cran.r-project.org/package=MendelianRandomization, 20 October 2020, date last accessed). We took all SNPs with BDNF GWAS P-value <1×10−7, and MAF > 0.01, and reduced these to a subset with pairwise r2 < 0.3 using the PLINK clumping tool (http://zzz.bwh.harvard.edu/plink/clump.shtml, 20 October 2020, date last accessed). For each SNP, the causal effect of BDNF levels on an outcome can be estimated as ratio of SNP effect size on the outcome to that on BDNF. A final causal estimate is a summary of the estimates based on individual SNPs, using the inverse-variance weighted method that accounts for LD between SNPs. MR-egg method was used to examine horizontal pleiotropy and to test for MR while accounting for potential horizontal pleiotropy, i.e. pathways for genetic variants affecting outcome other than through BDNF (Bowden ; Burgess and Thompson, 2017).

Co-localization

Using the coloc package, we tested whether each of the traits significantly associated with BDNF levels in the MR analyses shared similar associations with variants in the 1 Mb regions of the 7 loci associated with BDNF levels. Evidence of co-localization was evaluated using the posterior probability for hypothesis four (PP4), indicating that the BDNF levels and the other phenotypes shared the same variant association in the region. We set a prior probability of such association to 1×10−4.

Data availability

The genome-wide summary statistics that support the findings of this study will be made available via the database of Genotypes and Phenotypes (dbGaP) website (https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000930.v6.p1, 20 October 2020, date last accessed) upon publication. BDNF and genotype data are available from the corresponding authors Q.Y. and S.S. upon reasonable request.

Results

GWAS results

General characteristics of the study participants are provided in Supplementary Table 3. The genomic control factor λ= 0.86–1.03 across cohort indicates little or no inflation in the centre of P-value distributions that could be caused by stratification, cryptic relatedness or other violations of model assumptions (Supplementary Table 4). In addition, LD Score regression intercepts that measures inflation other than polygenicity ranged from 0.99 to 1.07. The ratio of intercepts over the mean chi-square statistic ranged from 0.98 to 1.004 indicating little inflation is due to polygenicity. The meta-analysis identified seven independent genome-wide significant (P-value < 5×10−8) loci (500 kb range, r2 < 0.3) associations (Table 1, Fig. 1 and Supplementary Fig. 1). Among them, six SNPs were located within genes and one was in intergenic region. The top-associated SNPs in the 7 loci were rs75945125 (P-value = 6.1×10−13, near BDNF), rs467369 (P-value = 2.6×10−11, BRD3), rs13084580 (P-value = 6.4×10−11, CSRNP1), rs3824987 (P-value = 2.2×10−10, KDELC2), rs71329093 (P-value = 1.2×10−9, RUNX1), rs2242882(P-value = 3.1×10−9, RUNX1) and rs1488831 (P-value = 3.6×10−9, BDNF-AS) (Table 1, Fig. 1A and B). According to the conditional and joint association analysis, the first 5 SNPs were independently associated with the trait. However, there was lack of association with rs6265 (P = 0.37) from looking up this SNP in the GWAS results. Hence, this SNP was not included in our conditional and joint analyses using GCTA-COJO. The results of all SNPs with P-value<10−7 and MAF>0.01 are presented in Supplementary Table 5.

Table 1

Genome-wide significant results from the meta-analyses of circulating BDNF levels

SNP	Chr:Pos	Gene	Distance	A1/A2	Freq1	Zscore	P	H2q
rs75945125	11:27778925	BDNF	trans, 56 kb	T/C	0.9533	−7.20	6.1 × 10−13	0.52%
rs467369	9:136905765	BRD3	cis, 0	T/C	0.5478	−6.67	2.6 × 10−11	0.45%
rs13084580	3:39188182	CSRNP1	cis, 0	T/C	0.1122	−6.55	6.4 × 10−11	0.40%
rs3824987	11:108346103	KDELC2	cis, 0	A/G	0.5928	−6.35	2.2 × 10−10	0.37%
rs71329093	21:36392564	RUNX1	cis, 0	A/G	0.9512	6.08	1.2 × 10−9	0.46%
rs2242882	21:36381815	RUNX1	cis, 0	T/C	0.1093	−5.93	3.1 × 10−9	0.33%
rs1488831	11:27637014	BDNF-AS	cis, 0	T/G	0.0648	5.90	3.6 × 10−9	0.31%

Seven independent variant with the lowest P-value in the fixed-effect effective-sample-size weighted Z-score meta-analysis is shown.

SNP, single-nucleotide polymorphism; Chr:Pos, chromosome and position; A1/A2, coded allele/non-coded allele; Freq1, frequency of coded allele; Z-score, Z-score from METAL; P, P-value; H2q, single SNP heritability.

Figure 1

Manhattan plots for GWAS meta-analyses of circulating BDNF. SNPs are plotted on the x-axis according to their position on each chromosome with -log10 P-value on the y-axis. The upper dashed horizontal line indicates the threshold for genome-wide significance, i.e. 5 × 10−8. Seven independent SNPs with the lowest P-value in the meta-analysis is coloured in orange.

We also utilized inverse-variance weighted meta-analysis to show the coded allele effect size of the lead variants. To harmonize the phenotypes across cohorts, we standardized the BDNF levels by its SD so that beta is the change in SD of circulating BDNF per increase of coded allele. No large variability was shown in beta estimates, and the test of heterogeneity of the betas across cohorts was not significant for all lead variants. Of note, the meta-analysed betas and SE from this analysis were very similar to the ones converted using the Z-score and MAF from sample size weighted meta-analyses (Supplementary Table 6, Fig. 2).

Figure 2

Forest plots in GWAS meta-analysis. Each plot denotes the association of the following SNPs with circulating BDNF levels: (A) rs75945125; (B) rs467369; (C) rs13084580; (D) rs3824987; (E) rs71329093; (F) rs2242882 and (G) rs1488831. Meta-analysis of beta were obtained by applying the inverse-variance weighted method. Column 1: Individual study name, summary indicates meta-analysis results. Column 2: The coded allele effect estimation, beta of individual study is the change in SD of BDNF per increase of coded allele. Column 3: The boxes show the effect estimates from the single studies, while the diamond shows the meta-analysis result. The horizontal lines through the boxes illustrate the length of the confidence interval. The width of the diamond serves the same purpose. The size of the box illustrates the sample size of each study. The vertical line is the line of no effect.

Our analyses on all FHS subjects showed that 30% (SE = 2.46%, P-value = 4×10−48) of variation in circulating BDNF concentrations can be attributable to genetic variance, suggesting that variability in circulating BDNF is moderately heritable. BDNF heritability for generation three subjects is 40% (SE = 4.37%, P-value = 3×10−26) and for the combined original cohort and offspring is 28% (SE = 5.04%, P-value = 2×10−9).

SNP heritability

The proportion of variance explained by the top seven SNPs derived from beta estimation and allele frequency ranged from 0.31% to 0.52%. Together they explained 2.84% of circulating BDNF variation (Table 1) in the data set they were derived in. We also estimated the SNP-based heritability of the lead variants from GWAS summary data accounting for the LD around the loci. Using HESS, we calculated the local SNP heritability of the six regions (60 kb range), containing the seven identified variants. The local SNP heritability estimation of the six region ranged from 0.31% to 0.57% and together they explained 2.62% variation in circulating BDNF (Supplementary Table 7).

Functional annotation

To functionally annotate our discoveries, we searched if the genome-wide association variants were also eQTL using data of PBMC, monocyte, brain and T cells from multiple sources and projects. We found that rs467369 is the eQTL of BRD3, rs3824987 is the eQTL of ACAT1, ATM and NPAT, rs13084580 is the eQTL of WDR48, TTC21A and SCN11A (Supplementary Table 8). We were also able to detect one protein quantitative trait loci. By associating the 7 SNPs to 228 ROSMAP proteins, we found rs75945125 is the protein quantitative trait loci of COX7B (beta =−0.50, P-value = 5.8×10−5, FDR = 0.09). Thus, functional annotation analysis allowed us to identify seven additional candidate genes (in addition to the six genes identified in our genome-wide meta-analysis, yielding a final set of 13 candidate genes; Supplementary Table 9). Furthermore, we related neuropathological findings to the expression of the 14 genes within 60 kb of GWAS lead variants in brains from the ROSMAP samples. Among the four phenotypes: neurofibrillary tangles, neuritic plaques, β-amyloid load and tau tangle density included in the analysis, the gene expression of BDNF was associated with less neurofibrillary tangles, using a Bonferroni threshold of 7×10−4 (Supplementary Table 10).

The module quantitative trait loci (modQTL) analysis fails to detect any significant result. Among the 47 modules of co-expressed genes from bulk brain, none was associated with the lead variants using an FDR threshold of 0.1. Similarly, by relating the seven lead SNPs to five ROSMAP phenotypes, no association was significant using an FDR threshold of 0.1.

Gene prioritization and pathway analysis

Based on FDR in DEPICT, no significant functional similarities among genes, expression enrichment across tissues, cell types or gene sets associated with circulating BDNF was identified.

Gene interaction networks

The list of 22 input genes located within 60 kb of the genome-wide significant SNPs is in Supplementary Table 11. A total of four networks were generated with scores 25, 16, 3 and 2, where the scores are -log10 (P-value) of the chance of observing such network with same number of randomly selected focus genes. The two most significant networks are presented in Figs 3 and 4. The most significant network (Fig. 3) contained 10 focus genes and 25 interconnecting molecules, with cell morphology, cellular function and maintenance, and embryonic development as top disease and functions. The second most significant network (Fig. 4) contained 7 focus genes and 28 interconnecting molecules, with cellular assembly and organization, cellular function and maintenance, cell morphology as top disease and functions. Canonical pathway analyses revealed that Ephrin-A Signaling, Granulocyte-macrophage colony-stimulating factor Signalling, Neurotrophin/TRK Signalling, pigment epithelium-derived factor Signalling, Acute Myeloid Leukemia Signalling, Neuropathic Pain Signalling in Dorsal Horn Neurons were among the most significant signal pathways.

Figure 3

The most significant interaction network from Ingenuity Pathway Analysis (Score = 25). Up-regulated molecules are in red-darker shade is more extreme.

Figure 4

The second most significant interaction network from Ingenuity Pathway Analysis (Score = 16). Up-regulated molecules are in red-darker shade is more extreme.

Two-sample MR analyses provided support for a potential causal role of BDNF on cardioembolism (CE; P = 0.007), white matter hyperintensity (P = 0.017) and high-density lipoprotein (HDL) cholesterol (P = 0.03) but not the other phenotypes considered (Table 2) at nominal significance level. However, after correcting for multiple testing involving seven domains (alpha = 0.05/7 = 0.007), only the MR test of CE remained significant. The associations of BDNF-related SNPs with CE are presented in Supplementary Table 12. The MR-egger analyses did not show nominal significant pleiotropy, as indicated by a non-significant intercept for the three traits (P > 0.05).

Table 2

Mendelian randomization analyses evaluating causal effects of BDNF on various phenotypes

Category	Outcome	Beta	SE	LL95CI	UL95CI	P-value	Reference
Cardiometabolic	Coronary heart disease	0.001	0.036	−0.070	0.071	0.988	Nikpay et al., (2015)
Cardiometabolic	Myocardial infarction	−0.048	0.038	−0.122	0.027	0.208	Nikpay et al., (2015)
Cardiometabolic	Cardioembolism	−0.376	0.141	−0.651	−0.101	0.007	Nikpay et al., (2015)
Cardiometabolic	Ischaemic stroke	−0.029	0.075	−0.176	0.117	0.695	Nikpay et al., (2015)
Cardiometabolic	Large vessel disease	0.161	0.166	−0.165	0.487	0.333	Nikpay et al., (2015)
Cardiometabolic	Small vessel disease	0.014	0.158	−0.295	0.324	0.928	Nikpay et al., (2015)
Cardiometabolic	Diabetes (BMI adjusted)	0.031	0.053	−0.074	0.135	0.565	Scott et al., (2017)
Cardiometabolic	Diabetes	0.016	0.045	−0.072	0.104	0.727	Scott et al., (2017)
Cognitive	General cognitive function	0.013	0.011	−0.009	0.036	0.24	Davies et al., (2018)
Cognitive	Alzheimer's Disease (IGAP 2019)	0.031	0.042	−0.052	0.113	0.468	Kunkle et al., (2019)
Cognitive	Alzheimer’s Disease (UKBB+IGAP 2013)	0.020	0.031	−0.041	0.081	0.523	Marioni et al., (2018)
Lipids	HDL	−0.063	0.029	−0.120	−0.006	0.030	Willer et al., (2013)
Lipids	LDL	0.016	0.032	−0.046	0.078	0.619	Willer et al., (2013)
Lipids	Total cholesterol	−0.003	0.031	−0.063	0.057	0.926	Willer et al., (2013)
Lipids	Triglycerides	0.010	0.029	−0.046	0.067	0.718	Willer et al., (2013)
Brain	hippocampus volume	−0.040	0.031	−0.101	0.022	0.204	Hibar et al., (2017)
Brain	White Matter Hyperintensities	−0.120	0.050	−0.214	−0.021	0.017	Verhaaren et al., (2015)
Lifestyle	Physical activity	0.014	0.022	−0.028	0.057	0.513	Doherty et al., (2018)
Lifestyle	Diet	0.021	0.025	−0.028	0.071	0.392	Cole et al., (2020)

The Mendelian randomization analyses were performed using two-sample inverse-variance method.

P-values that are bold have reached significance at P-value of <0.05.

Co-localization analysis

The posterior probabilities that CE and BDNF levels shared the associations with the same variant ranged from 0.15 to 0.43 across the seven loci, with the highest posterior probability obtained in the region of rs13084580. The highest posterior probability between BDNF and white matter hyperintensity was 0.43, at the regions of rs2242882 and rs71329093, respectively. There were two other regions with posterior probability >0.1 for white matter hyperintensity. The highest posterior probability between BDNF and HDL was 0.42 at the region of rs3824987, and no other posterior probability >0.1 for HDL was observed (Supplementary Table 13).

Sensitivity analysis

The GWAS meta-analysis among samples of serum BDNF only (excluding the Rotterdam Study cohort in which plasma levels were measured) resulted in similar findings (Supplementary Table 14). Among seven SNPs identified, three (rs467369, rs75945125 and rs3824987) were identified in the total sample. The other four SNPs (rs17791786, rs75967349, rs117571072 and rs2834730) were located within the same genes as in the main GWAS (RUNX1 and BDNF-AS) or in the eQTL analysis (TTC21A).

Discussion

In the current study, utilizing a sample from four cohort studies including a total of 11 785 individuals, mostly Caucasian, we have revealed previously unknown genetic architecture of circulating BDNF levels.

First, the current study shows that BDNF levels in human serum are heritable. The heritability of circulatory BDNF levels in our sample was moderate (∼30%), yet this estimate is high considering the strong environmental influences on blood BDNF concentrations such as exercise (Coelho ; Weinstein ) as well as other clinical and lifestyle factors (Weinstein ). Our heritability estimate is somewhat smaller compared to a previous estimate of 48% in a different study (Terracciano ). Yet, it should be noted that the current estimate of ∼30% may be more valid, as it is based on age- and sex-adjusted analysis utilizing extended pedigree rather than on unadjusted sib-pairs correlations in the previous study.

Our meta-analysis GWAS revealed seven novel SNPs significantly associated with serum BDNF levels. Three were located on chromosome 11, of them one maps near the BDNF gene, one in the BDNF-AS which transcribes a non-coding RNA and another is located within the KDELC2 gene and is protein coding. The other SNPs were scattered across chromosomes 3, 9 and 21 and were all mapped within protein-coding genes. While six out of the seven SNPs were located within genes, none was located within the BDNF gene. Rather, the SNP most significantly associated with circulating BDNF levels was rs75945125, which maps on chromosome 11 near the BDNF gene. Interestingly, the well-studied BNDF functional variant VAL66MET (rs6265) in the BDNF gene was unrelated to circulating BNDF levels in our study, a finding supported by some studies (Jiang ; Terracciano ) but not by others (Jin ; Kaess et al., 2015). These inconsistencies may be attributed to the complexity of the VAL66MET expression, which may be influenced by various factors including effect modification by age, sex, ethnicity and environment, as well as by the genetic model used for analysis, and gene–gene interactions [reviewed in Tsai (2018)]. Of note, a previous GWAS of circulating BDNF found some evidence of associations in polymorphism within the BDNF gene as well as near and within the NTRK3 gene which our study did not replicate (Terracciano ). Yet, it is important to note that the observed SNPs in this previous study did not reach statistical significance. In addition, the different findings in the previous study can be attributed to its smaller sample (N = 2054) and the inclusion of individuals from a genetically isolated population from Sardinia, Italy.

We additionally identified two SNPs within the RUNX1 gene that were linked with peripheral BDNF levels. The protein encoded by this gene represents the alpha subunit of core-binding factor and is thought to be involved in the development of normal haematopoiesis. The haematopoietic system is implicated in Alzheimer's disease (Ray ) and has been suggested as a promising target for Alzheimer's disease treatment (Lampron ). Dysfunctions within the haematopoietic system are thought to contribute to Alzheimer's disease aetiology through regulation of microglia homeostasis, a prominent feature of Alzheimer's disease (Hansen ). Indeed, evidence shows that the RUNX1 gene tightly regulates the development as well as the homeostasis of microglia cells (Kierdorf and Prinz, 2013). This role may be linked with a possible involvement in regulation of BDNF levels as evidence shows that BDNF can act as a direct chemotactic factor of bone-marrow-derived cells to promote recruitment of myeloid cells in cites of vascular injury (Kermani ), and in brain, it induces proliferation of microglial cells (Gomes ). On the other hand, BDNF is released from activated glia, and such ‘microglia BDNF’ is an important regulator of synaptic plasticity and function (Parkhurst ).

An eQTL analysis of the seven top SNPs reported additional genes potentially implicated in regulation of BDNF serum levels. Among them, the ACAT1 gene encodes a mitochondrially localized enzyme and is critical to maintain cholesterol homeostasis by converting free cholesterol to cholesteryl ester. Dysregulation of cholesterol pathways has been linked with increased Alzheimer's disease risk (Di Paolo and Kim, 2011), and indeed, findings from cell cultures and mouse models show that blocking the ACAT1 activity produces beneficial effects on Alzheimer's disease, such as decreasing amyloid-beta production (Shibuya ) as well as clearance of amyloid-beta from microglia (Shibuya ). In addition, a search using http://www.humanmine.org, 20 October 2020, date last accessed found that the protein product of both BDNF and BRD3 interact with β-amyloid (Aβ) (Oláh ).

Another gene identified in our eQTL analysis is the ATM gene. Although it was mostly linked with cancer risk in previous literature, it has been recently associated with coronary artery disease (Ding ) and type-2 diabetes (Ding ), both are major risk factors for cognitive impairment. The expression of another gene, SCN11A, has also been identified in Alzheimer’s disease transcriptome studies (Magistri ; Yang ) and has been linked with transmission of nerve impulse (Magistri ). Other genetic polymorphisms have been associated with serum BDNF concentrations in our GWAS meta-analysis (SNPs in BDNF-AS, BRD3, CSRNP1 and KDELC2), eQTL (SNPs in WDR48, TTC21A) and protein quantitative trait loci (an SNP in COX7B) analyses. Overall, future research is warranted to investigate the implications of these genes in brain aging, with particular interest to the BRD3 and CSRNP1 genes which in our expression study were related to cognitive decline.

MR conducted in our study pointed to a possible causal relationship of serum BDNF levels with CE. The link between BDNF and CE is supported by a growing body of evidence suggesting a potential involvement of BDNF in heart diseases (Ejiri ). In mice models, BDNF levels in periphery and brain have been shown to increase after induced myocardial infarction, yet genetic disruption of neuronal BDNF expression inhibited the increase of plasma BDNF after myocardial infarction and led to exacerbation of cardiac dysfunction (Okada ). In addition to experimental evidence, epidemiological studies also show that circulating BDNF levels are related to multiple risk factors for cardiovascular dysfunction (Golden ). Additionally, in previous research from the FHS cohorts, diminished circulating BDNF levels were observed in individuals who had a history of atrial fibrillation (Weinstein ), and in a separate study were prospectively associated with increased risk of cardiovascular disease (Kaess et al., 2015). While in the current analyses, there was no causal association between circulating BDNF and cardiovascular disease in general, our findings strongly suggest that BDNF plays a role in plaque instability. These findings are in line with accumulating evidence, showing that BDNF is implicated in angiogenesis and the maintenance of vascular integrity in recent reports (Donovan ). These angiogenic effects are attributed, at least partly, to its role as a direct chemotactic factor of bone-marrow-derived cells discussed above (Kermani and Hempstead, 2007). Further evidence comes from a former study demonstrating increased BDNF levels in coronary circulation of patients with unstable angina compared to those who suffered from stable angina, thus stressing the possible role of BDNF in plaque stability in line with results from animal models (Ejiri ). In the context of brain aging, plaque instability may serve as a possible pathophysiological mechanism. For example, evidence suggests that spontaneous cerebral emboli may represent a potentially treatable target to slow and possibly prevent dementia, both Alzheimer’s disease and vascular dementia (Purandare and Burns, 2009). Moreover, it has been shown that spontaneous cerebral emboli predict more rapid progression of dementia over 2 years in both Alzheimer’s disease and vascular dementia (Purandare ). Recent findings additionally suggest that measurement of carotid atherosclerotic plaque instability may be more clinically relevant to assess one’s cognitive function than a measurement of the degree of stenosis (Dempsey ). While in the context of these findings it is expected that BDNF serum levels may have shared genetic basis with cardioembolic stroke, previous outcomes included all ischaemic stroke subtypes (Malik ), and thus may lack the power to detect such associations.

Of note, the lack of causal association of BDNF levels with Alzheimer’s disease and general cognition in our MR may simply stress the complexity, multi-aetiological nature of the Alzheimer’s disease phenotype (Scheltens ). Hence, it is possible that these associations were not observed because complex interrelationships such as between gene (e.g. Val66Met) and environment (e.g. physical activity) were not considered (Brown ).

Circulating BDNF levels were associated with four significant networks with the lead motifs of cellular morphology and function. The most significant networks include additional genes and proteins that have been implicated in brain aging and Alzheimer’s disease, such as the Synapsin complex (Evergren ), the Metabotropic Glutamate Receptors (MGLURs) (Ribeiro ) and the extracellular signal-regulated kinase (ERK1/2) (Sun and Nan, 2017). The second network also included the ELAV which has been implicated in Alzheimer’s disease (Amadio ) and TRIM25, an E3 ubiquitin ligase, which is thought to be associated with acute and chronic neuroinflammation (Torre ). Furthermore, canonical networks highlighted several potential interactions between serum BDNF levels and molecular signalling. Of note, the Ephrin-A signalling has been implicated in brain development, synapse formation and plasticity (Lai and Ip, 2009) and controls brain size by regulating apoptosis of neural progenitor cells (Depaepe ). Similarly, the granulocyte-macrophage colony-stimulating factor signalling may have multiple neuroprotective effects in Alzheimer’s disease (Kiyota ) and brain injury (Shultz ), and is thought to play a major role in structural plasticity relevant to memory and learning (Krieger ). Lastly, the pigment epithelium-derived factor Signalling has been associated with cardiometabolic disorders (Yamagishi and Matsui, 2014) and development of atherosclerosis (Ma ).

This strength of our study is the substantial contribution it makes to our understanding of the genetic determinants and biological implications of variation in circulating BDNF levels, an intriguing molecule that affects many aspects of physiology ranging from memory, mood and cardiac function to smoking habits; BDNF may also be a biological mechanism through which physical and social activities, diet and other lifestyle factors influence disease risk. Limitations include the predominantly European ancestry of the persons studied and the restriction to less common and common variants. Further large-scale studies are warranted to identify rare variants at the loci we have identified, as well as epigenetic and gene–gene (G × G) or gene–environment (G × E) interactions, which may help explaining a greater proportion of the observed heritability. Lastly, there is a possibility of inflation in our MR study due to subtle population stratifications, assortive mating and dynastic effects (Hartwig ; Kong ), which may be overcome by utilizing data from large, family-based studies (Brumpton ).

Supplementary material

Supplementary material is available at Brain Communications online.

Click here for additional data file.