The brain-derived neurotrophic factor Val66Met polymorphism increases segregation of structural correlation networks in healthy adult brains

Participants

The Local Ethics Committee of the University of Occupational and Environmental Health, Kitakyushu, Japan (UOEH) granted Ethical approval to carry out the study (Ethical Application Ref: 第セH25-13号) A copy of the IRB approval documentation and the blank consent form are provided as supplemental files. All participants were informed about the purpose of the study, and written informed consents were obtained from all subjects via the forms.

Forty-nine subjects of Japanese origin were included in the study (mean age 41.1 years, standard deviation = 11.3 years; range 20–65 years; male:female = 36:13; see Table 1).

Most of the participants were recruited from the psychiatric or radiology departments of UOEH. None of the participants reported any neurological and/or psychopathological disorder (e.g., depression, attention deficit/hyperactivity disorder, and/or schizophrenia) in a simple questionnaire enquiring about their lifetime history of such diseases. All of the participants were right handed. We also asked for the participants’ number of years of education; this variable might help to better compare the results of future studies investigating the effect of the BDNF Val66Met polymorphism on the structure of the brain. However, as nearly all participants were medical doctors, no differences in education were observed. Because of the recruitment procedure used (medical doctors, characterized by a higher percentage of male doctors), males were over-represented in the investigated sample. All participants were of East Asian ethnicity.

After undergoing MRI scans, two radiologists with expertise in neuroradiology evaluated the images and confirmed that there were no abnormalities in the participants’ brain. Then each participant provided a blood sample for genotyping the BDNF Val66Met polymorphism.

Genetic analysis

DNA was extracted from the blood samples of each of the 49 participants according to standard laboratory protocols. DNA was isolated from peripheral blood mononuclear cells using the QIAamp (R) DNA Mini-Kit (QIAGEN, Tokyo, Japan). Genotyping was carried out with a polymerase chain reaction single nucleotide polymorphism (SNP) genotyping system using a BigDye Terminator v3.1 Cycle Sequencing Kit (Life Technologies Japan, Tokyo, Japan). The DNA was read using a BMG Applied Biosystem 3730xI DNA Analyzer (Life Technologies Japan, Tokyo, Japan). We used a forward primer (ATGAAGGCTGCCCCCATGAAA) and a reverse primer (TGACTACTGAGCATCACCCTG) for the BDNF Val66Met polymorphism. The participants were either homozygous for the Val allele (Val/Val genotype), heterozygous (Val/Met genotype), or homozygous for the Met allele (Met/Met genotype).

The 49 study subjects were divided into groups based on their BDNF genotype. Because the Met/Met genotype is less prevalent in Japan (about 15.9% (Shimizu, Hashimoto & Iyo, 2004)), participants were grouped according to the occurrence of the Met allele, which resulted in two independent groups: Val homozygotes (Val/Val genotype) and Met carriers (Val/Met and Met/Met genotypes).

Analyses of demographic and genotyping data

To assess the demographics and characteristics of participants, a nonparametric rank test (Wilcoxon–Mann–Whitney test) was applied to compare the differences in age and years of education among the subjects that were either Met carriers or Val homozygotes. Fisher’s exact test was used to evaluate the difference in gender among the groups. To assess the existence of any sampling bias within our subjects, we confirmed that the genotype distributions did not deviate from the expected frequency under Hardy–Weinberg equilibrium and tested the expected Japanese genotype/allele frequencies (Shimizu, Hashimoto & Iyo, 2004) of the BDNF gene polymorphism using chi-square tests.

These statistical analyses were performed using R software (version 3.4.1, R Foundation for Statistical Computing, Vienna, Austria). The level of significance was set to p < 0.05 for all demographic experiments.

MRI acquisition

MRI scans were performed on a 3.0-Tesla scanner (Signa EXCITE 3T; GE Healthcare, Milwaukee, WI, USA) using a dedicated eight-channel phased-array coil (USA Instruments Aurora, OH, USA). Three-dimensional fast spoiled gradient recalled acquisition with steady state (3D-FSPGR) images of the whole head were obtained with the following parameters:10/4.1/700 (repetition time ms/echo time ms/inversion time), a flip angle of 10, a 24- cm field of view and 1.2-mm-thick sections with 0.47 × 0.47 × 0.6 mm³ resolution. In addition, diffusion tensor images (DTIs) were also obtained assuming other hypothesis study. Although these DTIs were not used for this study’s analysis, the images were used to check whether the subject brains had any abnormalities.

Regional cortical thickness measurement

We used the FreeSurfer software package (version 6.0, https://surfer.nmr.mgh.harvard.edu/) for cortical surface reconstruction and cortical thickness estimation. Because detailed procedures for using FreeSurfer have previously been by authors such as (Dale, Fischl & Sereno, 1999; Fischl & Dale, 2000; Desikan et al., 2006; Balardin et al., 2015), we describe them here only briefly. The FreeSurfer processing stream consists of the co-registration of the subject T1 image to an atlas for bias correction and intensity normalization, and for brain extraction/removal of non-brain tissue. The cortical boundaries between the gray and white matter, and between the gray matter and cerebrospinal fluid, were tracked, tessellated and smoothed to produce a surface mesh. Topology correction and surface deformation were applied and cortical thickness was estimated by considering the closest distance from the gray/white matter boundary to the gray matter/cerebrospinal fluid boundary at each vertex on the tessellated surface. The cortical surface of each hemisphere was then parcellated based on gyral and sulcal landmarks according to the Desikan-Killiany atlas (Desikan et al., 2006). After FreeSurfer processing, 34 separated cerebral regions for each hemisphere (a total of 68 regions) were identified. The 68 brain regions’ names and the abbreviations for these names are summarized in Table S1.

Because a research member who was an expert in using FreeSurfer processing visually checked the final FreeSurfer outputs for accuracy or processing failures in a systematic manner, and all of the subjects’ outputs were confirmed to be correctly labeled and measured. Thus, we obtained 49 regional thickness values for each of the 68 brain regions.

Prior to the construction of structural correlation networks, FreeSurfer cortical thicknesses were corrected using the following steps: first, in order to remove the effects of age, gender, and an age-gender interaction, linear regression model analyses was conducted at every region with the best model based on the Akaike information criterion (AIC). Next, the residuals of the linear models were substituted for the raw cortical thickness values.

We examined the statistical differences in the cortical thickness of each region between the two groups in a Analysis S1. The results are summarized in Table S2.

Network construction

For network construction, we first defined anatomical connections as statistical associations in cortical thickness between the brain regions in order to characterize human brain networks. Such a morphometry-based connection concept has been described in many previous studies (Worsley et al., 2005; Lerch et al., 2006; He, Chen & Evans, 2007; Váša et al., 2018).

We evaluated statistical similarities in cortical thickness between all region pairs using Spearman correlation coefficients (r) across subjects in each group. The correlation coefficient r ranges in value from 1 (positive correlation) to −1 (negative correlation), and the closer to 1 or −1 the more significant the correlation. However, to avoid the complication of statistical feature descriptions in subsequent graph theoretical analyses, we converted r values into absolute values. (As a supplementary analysis, we conducted a negative-edge-removed version analysis to check if the results change when edges with negative weight were removed. Please refer to Analysis S1). In addition, we conducted edge selection based on Random Matrix Theory (RMT) as inspired by Deng et al. (2012). RMT was first proposed by Winger and Dyson as a powerful method to identify and model phase transitions associated with disorder and noise in statistical physics and material science (Wigner, 1967). The RMT method is able to automatically identify a threshold to construct a mathematically meaningful network from a correlation matrix. We used a RMT-based approach to conduct threshold scanning. Correlation coefficients greater than the threshold were retained, and the remaining coefficients were set to zero.

Thus, two inter-regional correlation matrices (N × N, where N is the number of brain regions, here N = 68) consisting of correlations between every pair of brain regions were acquired using the 16 subjects’ brains included in the Val homozygote group and 33 subjects’ brains included in the Met carrier group respectively.

To perform graph theoretical analyses, we used the anatomical inter-regional correlation matrices obtained above as the adjacency matrices of undirected weighted graphs. Nodes and edges represented brain regions and connections between the regions, respectively. The edges weight w_ij (an element of the adjacency matrix) corresponds to —r— between brain regions i and j; however, they were set as 1/—r— when calculating shortest path lengths.

Network comparison

For network comparison, we first visualized the comparisons between pairs of networks using heat maps and circular graphs. The heat maps were made with R-package ggplot2 (version 3.2.1; https://ggplot2.tidyverse.org/). The circular graphs were made using the Python software program (version 3.6.8; http://www.python.org) and the MNE-Python (v0.19) package (http://www.mne.tools/stable/index.html) (Gramfort et al., 2013).

Next, we compared the networks using network analysis. According to our hypothesis on the BDNF SNP’s effects on integration, segregation, and modularity, we examined the differences in average shortest path length (avgSPL), average local clustering coefficient (avgLCC), average local efficiency (avgLEFF), small-worldness (SW), and modularity (Q) inspired by previous studies (Latora & Marchiori, 2001; Rubinov & Sporns, 2010; Rubinov & Sporns, 2011). avgSPL is a representative measure for integration, and avgLCC and avgLEFF are well used measures of segregation. SW is a balance between avgSPL and avgLCC, and Q is a representative measure for the modular architecture of a network. We computed these network measures using R software version 3.3.1 (https://www.R-project.org) and the R-package igraph version 1.0.1 (http://www.igraph.org).

avgSPL is the average shortest path length among all reachable node pairs and was calculated using the function shortest.path in the igraph package.

avgLCC indicates the average strength of local-scale connectivity and is defined as $\frac{1}{N}{\sum }_{i=1}^N{C}_i^w$, where $C_i^w$ is a weighted version of the nodal clustering coefficient defined as $C_i^w=\frac{1}{s_i\left(k_i-1\right)}{\sum }_{j,h\in nn\left(i\right),i\ne j\ne h}\frac{w_{ij}+w_{ih}}{2}sign\left(w_{ij}{w}_{ih}{w}_{jh}\right)$. Here $nn\left(i\right)$ is the set of the nearest neighbors of node i. k_i is the degree of node i; $k_i={\sum }_{j=1}^{N}sign\left(w_{ij}\right)$, and $sign\left(x\right)$ is the sign function. $C_i^w$ was computed using the function transitivity in the igraph package.

avgLEFF indicates average information efficiency among nearest neighbor (i.e., at local level), and in is defined as $\frac{1}{N}.{\sum }_{i\in N}{E}_{loc}^w\left(i\right)$, where $E_{loc}^w\left(i\right)$ is the local efficiency of node i defined as $E_{loc}^w\left(i\right)=\frac{{\sum }_{j,h\in nn\left(i\right),j\ne h}{\left[d_{jh}^w\right]}^{-1}}{k_i\left(k_i-1\right)}$, where $d_{ij}^w$ is the shortest path length between i and j in a weighted network.

SW (Humphries, Gurney & Prescott, 2006) was proposed as a measure of the small-world property in real-world networked systems (Watts & Strogatz, 1998). The small-world property indicates that all node pairs in a network were reachable by a short distance (i.e., the distance expected from random networks), even though the network is divided into highly interconnected clusters (i.e., the network is far from a random network). Specifically, SW was calculated based on the avgSPL and avgLCC in the actual and randomized networks: (avgLCC_act/avgLCC _rand)/(avgSPL_act/avgSPL_rand), where X_act represents a network measure X (i.e., avgLCC or avgSPL) in the actual networks, and X_rand is the average X obtained from 100 randomized networks. For an actual network, the randomized networks were generated by permutation of the edge weights in the original network.

The modularity of networks was measured using the Q-value, an excellent guide to whether a particular division of a network into communities is strong or weak division (Newman, 2004a; Newman, 2004b). We considered a weighted version of the Q-value (Fortunato, 2010). The Q value is defined as the fraction of edge weights that lie within, rather than between, modules relative to that expected by chance as follows: ${\displaystyle Q=\frac{1}{2W}\sum _{ij}\left(w_{ij}-\frac{s_i{s}_j}{2W}\right)\delta \left(c_i,c_j\right),}$ where W is the sum of the weights of all edges, and $\delta \left(c_i,c_j\right)=1$ if nodes i and j belong to the same module but 0 otherwise. A network with a higher Q indicates a higher modular structure. In the present study, an algorithm based on simulated annealing (Guimerà & Nunes Amaral, 2005) was used to find the maximum Q in order to avoid the resolution limit problem in community (or module) detection (Fortunato & Barthélemy, 2007; Fortunato, 2010) wherever possible. The maximum Q was defined as the network modularity of brain networks. We used the netcarto function in the R-package rnetcarto (version 0.2.4) to compute network modularity Q.

We did not consider regional network characteristics (i.e. local network measures) because we could not find any report on the BDNF SNP’s effects on specific brain regions. To show statistically significant differences in the network measures between Val homozygotes and Met carriers, we calculated type I error rates using the permutation method under the null hypothesis. Details of the statistical method are described in the next section.

Network measure comparison with permutation method

To statistically test the difference in the network measures between the Met carriers’ network and the Val homozygotes’ network, we considered a permutation test approach for weighted correlation networks (Bassett et al., 2008; Gill, Datta & Datta, 2010). This approach randomly permutes the regional cortical thickness data sets from the subjects, and contrasts two networks built by partitioning the 49 subjects into two groups. We calculated the network measures from the permuted networks and obtained the approximate p-value by computing the probability that the simulated differences are greater than or equal to the observed difference. That is, ${\displaystyle p\left(DiffN{M}_{simulated\left(i\right)}>DiffN{M}_{observed}\right)=\frac{1}{tpc}\sum _{i=1}^{tpc}I\left(DiffN{M}_{simulated\left(i\right)}>DiffN{M}_{observed}\right)}$ ${\displaystyle p\left(DiffN{M}_{simulated\left(i\right)}<DiffN{M}_{observed}\right)=\frac{1}{tpc}\sum _{i=1}^{tpc}I\left(DiffN{M}_{simulated\left(i\right)}<DiffN{M}_{observed}\right)}$ where p(x) is the p-value under the condition x. tpc is the total permutation count (in this study, tpc=1000). DiffNM_simulated(i) is the difference between the network measures calculated from the permutated Met carriers’ network and permutated Val homozygotes’ network obtained at i th trial, and DiffNM_observed is the difference in the network measure between the observed Val homozygote and Met carrier networks. The p-values of every network measure were computed simultaneously using the same set of random permutations.

Under the null hypothesis the correlation pattern of each region is the same: the hypothesis was tested based on this permutation scheme being confirmed. That is, between the pseudo-two groups generated by permutation simulation, we calculated the network measure differences, and we examine the frequency that each simulated difference is greater (or less) than the observed difference. If the frequency is less than the predetermined significant level, then we reject the null hypothesis that there is no significant difference between the two groups, and we accept the alternative hypothesis that there is a significant difference between them. A one-sided nonparametric permutation test was used to test for global differences of the network measures between the two groups as specified by the hypotheses. The number of the repetition time for the permutation tests was set to 1,000 based on the rules of thumb and the load on the computer. We controlled the familywise error rate with the Bonferroni correction; the critical value (alpha) for an individual test by dividing the familywise error rate (0.05) by the number of tests. Since we conducted five statistical tests respecting global network measures in this study, the critical value for an individual test was 0.05/5 = 0.01, and we only consider individual tests with p-value < 0.01 to be significant. Statistical analyses were performed by using R software (version 3.6.1, R Foundation for Statistical Computing, Vienna, Austria).

Demographic and genotyping data

Table 1 displays the genotype frequency, gender, age, and years of education in each of the two groups. Age, years of education, and gender were not significantly different in the two genotype groups (p-values = 0.73, 0.996, and 0.508 respectively).

Furthermore, the genotype distribution of the Val66Met polymorphism in our healthy Japanese cohort comprised 16 Val/Val homozygotes (32.7%), 28 Val/Met heterozygotes (57.1%), and 5 Met/Met homozygous subjects (10.2%). The distribution of genotypes did not differ from the frequencies expected from Hardy-Weinberg equilibrium (χ² = 2.03, df = 1, p-value = 0.363). The allele and genotype frequencies are consistent with those reported previously for Japanese individuals (Shimizu, Hashimoto & Iyo, 2004) (allele: χ² = 0.109, df = 1, p-value = 0.741; genotype: χ² = 1.478, df = 2, p-value = 0.478).

Network comparison

First, to understand the overviews of the structural correlation networks of both groups, we constructed heat maps (Fig. 2). When the heat maps of the two groups were compared, a grid pattern was more clearly perceived in the Met carrier heat map than in the Val homozygote heat map. In a heat map of the brain network, the grids represent sub-divided brain regions (i.e., modules). From these grid patterns, we observed that the Met carrier brain network had clearer module architecture than the Val homozygote brain network. When we focused on the colors in each heat map, the Met carriers’ map had more red areas whereas the Val homozygote heat map had more green areas. This observation suggested that the Val homozygote networks consisted of middle connections whereas the Met carrier networks had more connections with high strength.

Next, we constructed circular graphs using the 100 strongest connections (Fig. 3). In both groups, each brain region had a strong correlation with the counterpart of the region which is represented by the white lines straddling the midline of the each circular graph. This finding seemed to reflect the anatomical left-right symmetry of the cerebrum. In addition, the Met carrier network had fewer brain regions related to the 100 strongest connections than the Val homozygote network. When we focused on the color of the edges, there is no apparent difference between the Val homozygote graph and Met carrier graph.

The results of the statistical comparisons of network measures related to integration, segregation (i.e., avgSPL, avgLCC, avgLEFF, and SW), and modular architecture (i.e. Q) are given in Table 2. The avgLEFF of the Met carriers’ network was larger than that of the Val homozygotes’ network, and the statistical difference between the avgLEFF of the two networks was significant (p-value < 0.01). Although the Met carriers network’s avgLCC and SW were larger, and avgSPL and Q were smaller, than those of the Val homozygote network, there were no significant differences between the two groups in regard to these measures (p-value ≧ 0.01). The permutation distributions and permutation p-values for these network measures are shown in Fig. 4.