Partitioning heritability analyses unveil the genetic architecture of human brain multidimensional functional connectivity patterns.

Resting-state functional connectivity profiles have been increasingly shown to be important endophenotypes that are tightly linked to human cognitive functions and psychiatric diseases, yet the genetic architecture of this multidimensional trait is barely understood. Using a unique sample of 1,704 unrelated, young and healthy Chinese Han individuals, we revealed a significant heritability of functional connectivity patterns in the whole brain and several subnetworks. We further proposed a partitioned heritability analysis for multidimensional functional connectivity patterns, which revealed the common and unique enrichment patterns of the genetic contributions to brain connectivity patterns for several gene sets linked to brain functions, including the genes expressed preferentially in the central nervous system and those associated with intelligence, educational attainment, attention-deficit/hyperactivity disorder, and schizophrenia. These results for the first time reveal the genetic architecture of multidimensional brain connectivity patterns across different networks and advance our understanding of the complex relationship between gene sets, neural networks, and behaviors.

INTRODUCTION

A central goal of human genetics is to understand the genetic architecture of complex traits (Evans & Keller, 2018; Timpson, Greenwood, Soranzo, Lawson, & Richards, 2018). It is generally recognized that the brain may serve as an important endophenotype for connecting genes and behaviors (Gottesman & Gould, 2003; Khadka et al., 2013; Meyer‐Lindenberg & Weinberger, 2006). For example, the global efficiency of brain connectivity, which is a widely used graphic theoretical measurement (Achard & Bullmore, 2007), can track and predict cognitive performance (Bassett et al., 2009) and patterns of diseases (Fornito, Zalesky, & Breakspear, 2015). Meanwhile, the multidimensional connectivity patterns could provide rich information that serves as a unique brain fingerprint to reliably identify the same person across various scan sessions (Finn et al., 2015). The multidimensional connectivity patterns reflect brain maturation (Kaufmann et al., 2017) and can predict fundamental cognitive functions such as intelligence (Dubois, Galdi, Paul, & Adolphs, 2018; Finn et al., 2015), creativity (Beaty et al., 2018), and attention (Rosenberg et al., 2016). It could also reliably differentiate healthy individuals from those suffering psychiatric conditions, such as attention‐deficit/hyperactivity disorder (ADHD; Rosenberg et al., 2016) and schizophrenia (SCZ; Yahata et al., 2016).

The genetic basis of this important endophenotype, however, is largely unknown. In an early study, Glahn et al. (2010) used extended pedigree samples to demonstrate the importance of genetic factors to the default mode network's functional connectivity. Later, twin studies revealed significant additive genetic effects on resting‐state networks (Colclough et al., 2017; Fornito et al., 2011; Ge, Holmes, Buckner, Smoller, & Sabuncu, 2017; van den Heuvel et al., 2013). Nevertheless, these twin studies could not uncover the specific genetic factors contributing to the multidimensional traits of the brain.

To overcome these limitations, genome‐wide complex traits analysis (GCTA) tools have been developed to estimate the single nucleotide polymorphisms (SNP)‐based heritability of a certain trait in a large‐scale unrelated population (Yang, Lee, Goddard, & Visscher, 2011), providing a lower bound of the narrow‐sense heritability estimated in pedigree studies. Furthermore, by partitioning the genes into meaningful clusters or pathways, researchers can examine the genetic architecture of these traits (Chen et al., 2015; P. H. Lee et al., 2016; Yang, Manolio, et al., 2011). Other studies have used polygenic scoring analysis (Torkamani, Wineinger, & Topol, 2018) or the linkage disequilibrium (LD) score regression method (Finucane et al., 2015) to estimate the aggregated contribution of multiple SNPs to certain traits. Although most existing studies have applied these methods to scalar traits, a novel method has been recently proposed to estimate the heritability of multidimensional traits, such as brain shape, using unrelated subjects (Ge et al., 2016). To date, no study has investigated the genetic architecture of multidimensional traits by combining genome‐wide SNP heritability and genomic partitioning strategies.

The current study aimed to investigate the polygenetic architecture of multidimensional functional connectivity traits using a group of young, healthy and unrelated Han Chinese individuals (n = 1,704, average age = 20.5 years, ranging from 16 to 29 years, with 1,453 of them in the range of 18–22 years). This limited age range is quite suitable for studying the effects of genetic factors on the brain functional connectome, as existing studies have shown that the distinctiveness of the brain functional connectome increases rapidly from 8 to 18 years (Kaufmann et al., 2017). In addition, the Han population is genetically more homogenous than the more widely studied Western populations, providing stronger statistical power for detecting genetic effects. To further understand the genetic architecture of each functional brain network, we focused on five gene sets that have been found to affect brain functions and extended the approach used in Ge et al. (2016) based on a genomic partitioning strategy (Yang, Manolio, et al., 2011) to examine the patterns of genetic enrichment in each brain network.

MATERIALS AND METHODS

Participants

The participants were recruited from the Cognitive Neurogenetic Study of Chinese Young Adults (CNSCYA) Project, which recruited over 2,500 Chinese Han young adults from Beijing and Chongqing, with a focus on understanding the genetic, neural, and cognitive mechanisms of human cognition. In our study, we restricted our analyses to 1,822 subjects with both genomic and resting‐state imaging data. One hundred eighteen subjects were excluded for the following reasons: 62 were removed due to large head motions (a translation greater than 3 mm in any direction, or a rotation greater than 3°), 29 were excluded due to a lack of whole‐brain coverage in the resting‐state scan, 6 were excluded for poor genomic quality (participants showing missing genotyped SNPs >5%), 10 were excluded for missing gender or age information, and 11 were excluded due to close genomic relationships (genetic relatedness >0.05, see below). As a result, a total of 1,704 unrelated samples (652 males, aged 16–29 years) with high‐quality genomic data and resting‐state fMRI scans were included in our analysis. All subjects were healthy, free of any psychiatric conditions, and with normal or corrected‐to‐normal vision. Written informed consent was obtained from each participant after a full explanation of the study procedure. This study was approved by the Institutional Review Boards (IRBs) of Beijing Normal University and Southwest University, China.

Genotype quality control and imputation

A 4 ml venous blood sample was collected from each subject. Genomic DNA was extracted according to standard methods. Samples were genotyped on Affymetrix 6.0, Illumina OmniExpress, Illumina Zhonghua, or Illumina Omni2.5. Standard genome‐wide association quality control filters were applied to each data set individually using the Plink 1.9 (Chang et al., 2015; https://www.cog-genomics.org/plink2). We excluded SNPs with any of the following conditions: A minor allele frequency (MAF) of <5%, a per‐SNP missingness >5%, or a failing of the Hardy–Weinberg equilibrium test (p < 1 × 10−6). After that, six participants showing missing SNPs >5% were excluded from subsequent analyses (Table ). The cleaned genotype data were then recoded to the forward strand using functions and strand alignment data files produced and hosted by Will Rayner (https://www.well.ox.ac.uk/~wrayner/strand/index.html). Finally, we used the HRC‐1000G‐check‐bim.pl script to perform some final filtering and split data by chromosome (https://www.well.ox.ac.uk/~wrayner/tools/).

The preprocessed genotype data were then imputed against the 1,000 Genomes reference panel (phase 3 version 5, East Asian population) by using the Michigan Imputation Server (Das et al., 2016; https://imputationserver.sph.umich.edu/index.html). Previous studies have shown that for the Han Chinese population, 1000G phase 3 reference panel performs better than the Haplotype Reference Consortium (HRC) in terms of both the number of well‐imputed SNPs and imputation quality (Lin et al., 2018) and the R 2 between experimental and imputed genotypes was >95% for common SNPs (Genomes Project et al., 2015). We used Eagle (version 2.4) for phasing (Loh et al., 2016) and Minimac4 for imputation. The imputed SNPs from different data sets were then combined. SNPs meeting any of the following criteria were excluded: Imputation information score R 2 < 0.3 (Y. Li, Willer, Ding, Scheet, & Abecasis, 2010), Hardy–Weinberg p‐value < 1 × 10−6, MAF < 1%, or per‐SNP missingness >5%. This yielded ~6.6 million common autosomal chromosome SNPs after quality control.

To remove close relatives in the sample, we used GCTA functions (Yang, Lee, et al., 2011; https://cnsgenomics.com/software/gcta) to estimate the genetic relatedness for pairs of individuals in the combined data set using all of the ~6.6 million common SNPs. One sample from each pair of individuals with an estimated genetic relatedness >0.05 was removed (11 subjects were removed).

Furthermore, we performed principal component (PC) analysis (Price et al., 2006) in the combined data set using GCTA functions (Yang, Lee, et al., 2011) to check the ancestry of the participants used in the current study. Briefly, we used the overlapping SNPs (~3.3 million) of our current genome data and the 1000 genome project phase 3 data (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/), and the results show that all our participants are East Asians (Figure S1).

Image acquisition

Neuroimaging data were acquired with a 3.0 T Siemens MRI Trio scanner in the Brain Imaging Centers at Beijing Normal University and Southwest University. Anatomical MRI scans were acquired using a T1‐weighted, three‐dimensional, gradient‐echo pulse sequence. Parameters for this sequence were as follows: Repetition time/echo time/θ = 2,530 ms/3.39 ms/7°, field of view = 256 × 256 mm, matrix = 192 × 256, and slice thickness = 1.33 mm. A total of 128 and 144 sagittal slices were acquired to provide a high‐resolution structural image of the whole brain for the Beijing data scanned during phase one (2006–2008) and phase two (2013–2017), respectively. For Chongqing sample, repetition time/echo time/θ = 2600 msec/3.02msec/8°, field of view = 256 × 256 mm, matrix = 256 × 256, and slice thickness = 1.00 mm. A total of 176 sagittal slices were acquired.

For the resting‐state scan, participants laid supine on the scanner bed and were instructed to close their eyes and to not think about anything in particular. Foam pads were used to minimize head motion. Functional scanning used a gradient echo EPI sequence with PACE (prospective acquisition correction).

The following parameters were used: TR = 2,000 ms; TE = 30 ms; flip angle = 90°; FOV = 200 × 200 mm2; 64 × 64 matrix size with a resolution of 3.1 × 3.1 mm2. Thirty‐three 3.0 mm and 3.5 mm transverse slices were used for the Beijing data scanned during phase one (2006–2008) and phase two (2013–2017), respectively. The Chongqing sample used thirty‐two 3.0 mm transverse slices to cover the whole cerebrum and most of the cerebellum, TR= 2000 ms; TE = 30 ms; flip angle = 90°; FOV = 220 × 220 mm2; 64 ×64 matrix size with a resolution of 3.4 × 3.4 mm2. Two hundred forty brain volumes (time points) were acquired from the Beijing phase one sample, 200 volumes from the Beijing phase two sample, and 242 volumes from the Chongqing sample.

Image processing

Image preprocessing was performed using GRETNA tools (J. Wang, Wang, et al., 2015) and the AFNI toolbox (Cox, 1996). The first three volumes of each session were discarded by the scanner automatically to allow for signal equilibrium. In addition, the first 10 EPI volumes of the resting‐state scan were deleted to allow for signal equilibrium and to allow the participants to adapt to the scanning noise. The remaining images were slice‐time corrected, realigned, and registered to the standardized MNI space. Each fMRI volume was then segmented into separate tissue types (gray matter, white matter, and cerebrospinal fluid) using DARTEL (Ashburner, 2007). We then performed temporal detrending, nuisance regression, and bandpass filtering using AFNI tools. Based on the methods from a recent study (Lindquist, Geuter, Wager, & Caffo, 2019), we combined nuisance covariates, linear trends, and temporal filters (0.01–0.1 Hz) into a single regression model to avoid reinjection of the noise signals. The nuisance covariates included the average signal from the cerebrospinal fluid and white matter, the global signal, and 24 motion parameters (Friston, Williams, Howard, Frackowiak, & Robert, 1996). Finally, the data were spatially smoothed using a 4‐mm full‐width‐at‐half‐maximum (FWHM) Gaussian kernel.

Network construction

Functional connectivity was assessed using a widely used parcellation scheme consisting of 264 nodes across the whole brain (Power et al., 2011). We extracted the time series from each node by taking the mean signal in all voxels. Functional connectivity was calculated as the pairwise correlation between the time series from all nodes which resulted in a 264 × 264 connectivity matrix for each participant. We calculated the multidimensional functional connectivity patterns in the whole brain (all), as well as in the following 10 major subnetworks (Cole et al., 2013): Motor and somatosensory (somatomotor), cingulo‐opercular, auditory, default mode, visual, frontal–parietal, salience, subcortical, ventral attention, and dorsal attention networks.

Multidimensional functional connectivity

Pearson correlations between pairs of node time series were calculated and transformed using Fisher's Z‐transformation. The connectivity strength was then spatially standardized into Z scores to make the functional connectivity comparable across subjects (Ge et al., 2017). All the resulting values within each network were defined as the multidimensional functional connectivity.

Genome‐wide SNP‐based heritability analysis

SNP‐based heritability is defined as the proportion of phenotypic variance in a population that can be explained by the additive effects of a set of SNPs. We used methods recently developed by Ge et al. (2016) to estimate the genome‐wide heritability of the multidimensional functional connectivity patterns. Briefly, a given multidimensional trait Y, where Y is a matrix comprising M pair‐node time series correlations from N individuals, could be modeled as the following form: Y = G + C + E, where G, C and E are also N × M matrices, G represents the sum of additive genetic effects, C represents the shared environmental factors, and E represents measurement errors and other known sources of variance. Since our sample consisted of unrelated subjects, the shared environmental factors could be ignored, and the model can be redefined as Y = G + E. Where vec(G) ~ N (0, ∑A⊗K), vec(E) ~ N (0, ∑E⊗I), vec(.) is the matrix vectorization operator, ⊗ the Kronecker product of matrices. ∑A the genetic covariance matrix, ∑E the residuals covariance matrix, I an identity matrix, and K the genetic relationship matrix (GRM) for each pair of individuals estimated from genetic data, which was calculated by using GCTA functions (Yang, Lee, et al., 2011). Briefly, the genetic relationship between individuals a and b is defined as the mean correlation of SNP values between a and b over a number of SNPs:where denotes the number of copies of the reference allele (usually the minor allele) for the th SNP of the ath individual, and is the frequency of the reference allele.

The heritability of a set of multidimensional traits can be defined as:where ∑P is the phenotypic covariance matrix, and tr[.] is the trace operator for a matrix. Finally, the model is fitted using a moment‐matching estimator.

To adjust for the remaining subpopulation stratification, we performed the PC analysis using all of the ~6.6 million imputed data, and then included the top 10 principal components of ancestry, as well as age, gender, fMRI scanner, genotype array, and a head motion measure (mean framewise displacement) as covariates in the model.

Following Ge et al. (2016), the current study used both the parametric Wald test and the nonparametric permutation test to examine the statistical significance of genome‐wide SNP heritability. For the permutation test, we simultaneously shuffled the columns and rows of the GRM 10,000 times and recorded the heritability values estimated from the permuted data. Note that this approach is equivalent to randomly permuting the subjects in our study. Then, the estimated heritability from the original data was compared to the distribution of the permuted heritability values. Correction for multiple comparisons was performed with the false discovery rate (FDR; Benjamini & Hochberg, 1995).

Estimating enrichment of candidate gene sets

To better characterize the genetic architecture of functional connectivity, we further examined the enrichment (higher contribution) of heritability of certain gene sets in each brain network. We included five functional gene sets from previous studies: Genes preferentially expressed in the central nervous system (CNS; S. H. Lee et al., 2012; Raychaudhuri et al., 2010), SNPs associated with human intelligence (Savage et al., 2018), SNPs associated with educational attainment (J. J. Lee et al., 2018), SCZ‐associated SNPs (Z. Li et al., 2017), and ADHD‐associated SNPs (Demontis et al., 2018). We also included SNPs associated with Crohn's disease as a negative control set (Liu et al., 2015), which should show no association with brain connectivity patterns.

For the CNS genes, we defined the genic boundaries as 50 kb upstream and downstream from the 3′ and 5′ untranslated regions (UTRs) of each gene according to the UCSC hg19 assembly, which resulted in a set of ~1.3 million SNPs (~20% of the whole genome SNPs). We obtained the genome‐wide association studies (GWAS) summary results for the two psychiatric diseases from the Psychiatric Genomics Consortium (https://www.med.unc.edu/pgc/results-and-downloads/), the GWAS results for intelligence from the Complex Trait Genetics Lab (https://ctg.cncr.nl/), the GWAS results for educational attainment from the Social Science Genetic Association Consortium (https://www.thessgac.org/data), and the GWAS results for Crohn's disease from International Inflammatory Bowel Disease Genetics Consortium (https://www.ibdgenetics.org/downloads.html). All of the imputed SNPs in our current study were ranked based on p values from the GWAS results mentioned above, and we defined SNPs that passed a given threshold (here we selected the top 10%) as associated SNPs of the corresponding traits.

We divided our imputation data into two sets: One contained the trait‐associated SNPs and the other contained what were defined as control SNPs. Two genetic relationship matrices (GRMs) were then calculated using GCTA functions (Yang, Lee, et al., 2011). Previously, Yang et al. demonstrated that a joint analysis of multiple GRMs as random effects prevented inflated estimates of local SNP heritability (Yang, Manolio, et al., 2011). Thus, we used a joint analysis of the two genetic components to estimate the genome‐wide SNP heritability as the sum of (the heritability attributed to the trait‐associated SNPs) and (the heritability attributed to the set of unselected SNPs).

To examine gene enrichment of a set of genes for multidimensional traits, we partitioned the genetic covariance matrix ∑ into two matrices: ∑, the genetic covariance due to the candidate gene set, and ∑, the genetic covariance due to the unselected SNPs. Thus, the candidate‐set SNP heritability can be formulated as:

As in the genome‐wide SNP heritability analysis, age, gender, fMRI scanner, mean framewise displacement, genotype array, and top 10 principal components were included as covariates. Similar to Ge et al. (2016), we used a data transformation strategy to account for these covariates and a moment‐matching method to estimate the genetic covariances. Specifically, we estimated the standard error (SE) of the candidate SNP sets (SE[set]) using the jackknife resampling technique (Efron & Stein, 1981).

Enrichment was computed for each gene set as the ratio of the estimated to the expected , which was the genome‐wide heritability multiplied by the percentage of the SNPs in the given set among all SNPs.

The degree of enrichment was quantified by a Z‐score calculated as follows:

The p‐value for the Z‐score was calculated by assuming that its null distribution follows the standard normal distribution. Detailed method descriptions can be found in Gusev et al. (2014).

Correction for multiple comparisons was again performed with the FDR.

The method (using Z‐score to test for significance) is often used to test the significance of gene set enrichment, such as the widely used LD‐score regression (Finucane et al., 2015), and other studies (Tansey & Hill, 2018).

We did not use the nonparametric permutation test because it is computationally very intensive.

It is worth noting that previous studies often set heritability estimates between zero and one to exclude negative values; such a practice might lead to biased estimates, especially for moderate sample sizes (Wang, Yandell, & Rutledge, 1992). In order to obtain unbiased estimates of h 2, we allowed for negative values (i.e., did not set heritability estimates between zero and one), even though they were very rare.

Validation analysis

We performed several validation analyses to verify the robustness of our results. Specifically, we considered different parcellation schemes, different preprocessing parameters (i.e., global signal regression), and different numbers of selected candidate‐set SNPs.

Effects of parcellation scheme

To examine whether our main results could be replicated across different brain parcellation schemes, we performed genome‐wide SNP heritability analysis using a different brain functional parcellation scheme (i.e., Shen 268 nodes; Shen, Tokoglu, Papademetris, & Constable, 2013), which included the following 9 subnetworks: Medial frontal, frontoparietal, default mode, subcortical‐cerebellum, motor, visual 1, visual 2, and visual association, as well as executive control networks (the combination of medial frontal and frontoparietal networks). Only 1,533 subjects were eligible for this parcellation because the other subjects lacked whole‐brain coverage.

Effects of global signal regression

Global signal regression is still a controversial issue in resting‐state fMRI. Several studies have shown that a 24‐parameter head motion regression is not sufficient to remove motion confounds from functional connectivity estimates (Ciric et al., 2017; Parkes, Fulcher, Yucel, & Fornito, 2018), but global signal regression can improve the performance of this model. For this reason, we performed global signal regression in the main analysis. Nevertheless, other studies have suggested that global signal regression may alter network properties by introducing anti‐correlations (Murphy, Birn, Handwerker, Jones, & Bandettini, 2009; Weissenbacher et al., 2009). Recently, Murphy and Fox (2017) suggested that different preprocessing strategies are likely to provide complementary insights into functional brain organization, so we examined genome‐wide heritability without global signal regression in our validation analyses.

Effects of the number of SNPs in enrichment analysis

For the candidate‐sets analyses (i.e., SNPs associated with intelligence, educational attainment, SCZ, ADHD, and Crohn's disease), different numbers of SNPs were included depending on the threshold used. In the validation analyses, we also lowered the significance threshold to the top 20% and top 30% SNPs, yielding more SNPs in each candidate gene set.

RESULTS

Genome‐wide SNP heritability of multidimensional functional connectivity patterns

Our genome‐wide SNP‐based heritability analysis revealed that the whole‐brain (all) connectivity had significant heritability ( = 0.049, SE = 0.011, p Wald = 1.69 × 10−5, p permutation = 5.50 × 10−4, corrected). Four of the 10 subnetworks also showed significant heritability after FDR correction (Figure 1), including the visual network ( = 0.087, SE = 0.028, p Wald = 3.12 × 10−3, p permutation = 1.47 × 10−3, corrected), the default mode network ( = 0.065, SE = 0.015, p Wald = 4.23 × 10−5, p permutation = 5.50 × 10−4, corrected), the somatomotor network ( = 0.055, SE = 0.021, p Wald = 1.27 × 10−2, p permutation = 7.98 × 10−3, corrected), and the frontal–parietal network ( = 0.052, SE = 0.021, p Wald = 1.72 × 10−2, p permutation = 1.23 × 10−2, corrected).

FIGURE 1

Genome‐wide SNP heritability of multidimensional functional connectivity patterns across networks. Error bars represent the standard errors (SE) of the heritability estimates. Networks showing significant heritability after FDR correction are highlighted with asterisks (***p < .001, **p < .01, *p < .05)

Enrichment of candidate gene sets in the heritability of multidimensional functional connectivity patterns

The above results suggested that some of the functional networks in the current study were significantly heritable, making them appropriate targets for more specific genetic analyses. To further examine whether certain gene sets showed a relatively higher contribution (i.e., enrichment) to the heritability of multidimensional functional connectivity patterns, we selected five gene sets that have been linked to cognitive function, educational attachment, psychiatric diseases, and brain gene expression (Section 2), as well as SNPs associated with Crohn's disease, which served as a nonpsychiatric control gene set.

Existing methods have focused on the contribution of all available genetic variants to univariate or multidimensional traits (Ge et al., 2016; Yang et al., 2010). Here, we proposed a partitioned heritability analysis for multidimensional traits (Section 2). Using this novel approach, we discovered several interesting findings (Figure 2). In particular, we found genes that preferentially expressed in the central nervous system were enriched in almost all the networks that showed significant genome‐wide heritability, including the frontal–parietal network (2.71× enrichment, p = 1.97 × 10−91, corrected), visual network (2.17× enrichment, p = 1.46 × 10−25, corrected), default mode network (1.62× enrichment, p = 6.60 × 10−28, corrected), and the whole‐brain (all) network (1.57× enrichment, p = 2.02 × 10−48, corrected).

FIGURE 2

Enrichment of SNP‐heritability for multidimensional functional connectivity across candidate gene sets and networks. Each solid bar displays estimated enrichment folds, while the null hypothesis of 1.0 enrichment is shown by a dashed dark line. The p values indicate the significance of the difference from the expectation. Error bars represent the enrichment SE (= SE(set)/ % SNP). Networks with significantly enriched candidate set‐based heritability after FDR correction are marked with asterisks (***p < .001, **p < 0.01, *p < .05). We selected the top 10% SNPs of the ranked imputed genome data as trait‐associated variants. CNS, genes expressed in the central nervous system; IQ, SNPs associated with human intelligence; EA, SNPs associated with educational attainment; SCZ, SNPs associated with schizophrenia; ADHD, SNPs associated with attention‐deficit/hyperactivity disorder; Crohn, SNPs associated with Crohn's disease

Interestingly, we found that the candidate SNP sets associated with intelligence were enriched in all networks that showed significant genome‐wide SNP heritability. In particular, the intelligence‐associated genes were enriched in the visual network (4.68× enrichment, p = 5.43 × 10−93, corrected), whole‐brain (all) network (2.07× enrichment, p = 4.51 × 10−60, corrected), default mode network (2.29× enrichment, p = 2.44 × 10−43, corrected), the frontal–parietal network (1.50× enrichment, p = 2.68 × 10−4, corrected), and the somatomotor network (1.48× enrichment, p = 3.56 × 10−4, corrected). We also found that the SNP set associated with educational attainment showed enrichment in the visual network (3.55× enrichment, p = 5.54 × 10−54, corrected) and the frontal–parietal network (1.79× enrichment, p = 7.43 × 10−10, corrected).

For the two gene sets related to psychiatric diseases, SCZ‐associated SNPs showed enrichment in the whole brain (1.89× enrichment, p = 6.54 × 10−49, corrected), the frontal–parietal network (1.40× enrichment, p = 2.78 × 10−3, corrected), and the visual network (1.29× enrichment, p = 8.61 × 10−3, corrected). ADHD‐associated SNPs were also enriched in the frontal–parietal network (2.43× enrichment, p = 2.15 × 10−24, corrected) and the somatomotor network (2.42× enrichment, p = 1.77 × 10−26, corrected).

As expected, gene sets associated with Crohn's disease showed no enrichment in any of the above networks.

Validation results

Using the 268 parcellation scheme, we found that 5 of the 9 subnetworks and the whole‐brain (all) network were significantly heritable, whereas the default mode network showed marginally significant heritability ( = 0.062, SE = 0.037, p = .065, corrected; Figure 3). These results suggested that these brain networks, including the motor, visual, frontal–parietal, default mode, and whole‐brain networks showed reliable heritability and that the results were not affected by parcellation scheme. Without global signal regression, we also obtained very similar genome‐wide heritability results (Figure 4), except that the frontal–parietal network ( = 0.039, SE = 0.022, p = .107, corrected) and the somatomotor network ( = 0.033, SE = 0.023, p = .144, corrected) showed marginally significant heritability in this validation analysis. Finally, our results showed that different thresholds (and hence the different numbers of SNPs) generated similar enrichment patterns (Figure 5). For example, regardless of the threshold, the intelligence‐associated SNPs still showed enrichment in the default mode network, visual network, and the whole brain (all) network. Similarly, ADHD‐ and SCZ‐associated SNPs consistently showed enrichment in the frontal–parietal network for different threshold levels. SCZ‐associated SNPs were also robustly enriched in the whole brain network. Notably, with increasingly stricter threshold, we obtained higher folds of enrichment, suggesting the SNPs’ contributions to the phenotypes and their heritability were reliable.

FIGURE 3

Genome‐wide SNP heritability of multidimensional functional connectivity patterns across networks using the 268 nodes parcellation scheme. Global signal regression was applied. Error bars represent the SE of the heritability estimates. Networks showing significant heritability after FDR correction are highlighted with asterisks (***p < .001, **p < .01, *p < .05)

FIGURE 4

Effects of global signal regression on the SNP heritability of multidimensional functional connectivity patterns. Error bars represent the SE of the heritability estimates. withGSR, Global signal regression was applied in the preprocessing analysis. withoutGSR, Global signal regression was not applied. Networks showing significant heritability after FDR correction are highlighted with asterisks (***p < .001, **p < .01, *p < .05)

FIGURE 5

Effects of the number of selected associated SNPs on the enrichment pattern of multidimensional functional connectivity. The enrichment pattern overall was very similar across different numbers of SNPs included among the candidate sets. Each solid bar displays estimated enrichment folds, while the baseline (1.0 enrichment) is shown by a dashed dark line. Error bars represent the enrichment SE (enrichment SE = SE(set)/ % SNP). Reported are the results using three different cut‐off thresholds: Top 10%, top 20%, and top 30%. Networks with significantly enriched candidate‐set‐based heritability after FDR correction are highlighted with asterisks (***p < .001, **p < .01, *p < .05). See Figure 2 for definition of acronyms

Together, the above validation results suggested that our main results were robust and replicable.

DISCUSSION

To the best of our knowledge, this is the first study that examined SNP‐based heritability analysis and the genetic architecture of multidimensional resting‐state functional connectivity patterns across different networks in a large unrelated population using partitioned heritability methods (gene enrichment). Compared to univariate analysis, it has been suggested that joint multivariate analysis of multivariate traits could better capture the nature of data and increase statistical power (Ge et al., 2016; Schmitz, Cherny, & Fulker, 1998), which is a very helpful approach when the sample size is limited. In addition, partitioned heritability analysis revealed a significant shared genetic basis between certain functional networks and psychiatric diseases.

Using a very large sample size of an unrelated population from Biobank, Elliott et al. (2018) performed a comprehensive genetic association analysis of a number of brain imaging phenotypes. By examining each functional connectivity edge individually, they found that only 235 of 1,771 edges were significantly heritable; this proportion is much smaller than the proportion of networks showing significant heritability (5 of 11) in the current study. One factor contributing to the difference is that our joint analysis of multidimensional traits could increase statistical power (Schmitz et al., 1998). In addition, our sample was more homogenous in terms of age (young adults vs. middle‐aged to elderly adults) and genetic ancestry (Chinese Han vs. British), which could also improve statistical power of heritability analysis. Moreover, following Finucane et al. (2015), Elliott et al. (2018) partitioned the heritability into 24 functional categories such as coding, UTR, promoter, intron, enhancer, and so on, and revealed no enrichment for functional connectivity edges in any of these functional categories. In contrast, in the current study, we chose the candidate gene sets that have been associated with psychiatric diseases or cognitive functions, and obtained some meaningful gene set enrichment results. The different candidate gene sets used in the two studies are likely to contribute to the different enrichment results.

Our SNP‐based heritability estimate of multiple traits, which was a weighted average of the individual trait's heritability (Ge et al., 2016), was overall lower than that using twins from the Human Connectome Project (Colclough et al., 2017; where the average heritability of all connections was 15%). There are two potential reasons for this discrepancy. First, the genome‐wide SNP heritability in our study only captured additive genetic effects due to autosomal common SNPs, and hence ignored the effects of rare variants and sex chromosome variants and provided a lower‐bound estimate of the narrow‐sense heritability. Second, traditional twin and family designs may inflate the heritability estimation because they violate the shared common environment (Nolte et al., 2017).

Our partitioned heritability analysis revealed several important patterns of gene set enrichments that could help advance our understanding of brain‐behavioral relationships. First, we found that genes preferentially expressed in the CNS showed enrichment in almost all the networks showing significant SNP heritability. Previous studies have suggested that such CNS genes are associated with susceptibility to psychiatric diseases (Cross‐Disorder Group of the Psychiatric Genomics et al., 2013; S. H. Lee et al., 2012). Our results suggested that the CNS genes might contribute to psychiatric diseases by modulating the function of these networks, thus furthering our understanding of the biological mechanism of these genes.

Our study indicated that the set of SNPs associated with intelligence were enriched in all of the five networks showing significant genome‐wide heritability, including the default mode, frontal–parietal, visual, somatomotor, and whole‐brain networks. These networks are known to be involved in human intelligence (Deary, Penke, & Johnson, 2010; Ming et al., 2009; Sripada et al., 2019; Woolgar et al., 2010). In particular, the default mode network plays an important role in working memory, language, and intelligence (Greicius, Krasnow, Reiss, & Menon, 2003; Schultz & Cole, 2016; Smith, Mitchell, & Duncan, 2018). Furthermore, the posterior cingulate cortex, a key node in the default mode network, has been found to be a strong focus of cross‐network interactions (de Pasquale et al., 2012); thus, it may serve as a central and flexible network hub and play an important role in human cognition and related neuropsychological diseases (Leech & Sharp, 2014). Interestingly, although the set of SNPs associated with educational attainment also showed enrichment in the frontal–parietal network and the visual network, the overall enrichment was weaker compared with the set of SNPs associated with intelligence. This is consistent with the finding that educational attainment is a complex trait that is influenced by intelligence as well as various other factors such as personality, family income, school quality, and so on (Krapohl et al., 2014; Krapohl & Plomin, 2016).

We have observed both common and distinct patterns of enrichment for gene sets associated with the two different psychiatric disorders (ADHD and SCZ). In particular, both sets of SNPs showed enrichment in the frontal–parietal network. SCZ‐associated SNPs also showed enrichment in the visual network and the whole‐brain network, whereas only ADHD‐associated SNPs showed enrichment in the somatomotor network. These findings are consistent with the notion that common and distinct neural mechanisms exist among different psychiatric diseases (Khadka et al., 2013).

One interesting network associated with ADHD‐associated SNPs was the somatomotor network. Recent review and meta‐analysis articles note that a wide range of brain networks, including the frontal–parietal and somatomotor networks, are related to ADHD (Castellanos & Proal, 2012; Cortese et al., 2012). Furthermore, using transcranial magnetic stimulation, researchers have found abnormal inhibition in the motor network in ADHD patients (Gilbert, Isaacs, Augusta, MacNeil, & Mostofsky, 2011). The somatomotor network may act as a compensatory mechanism in patients with ADHD, who usually exhibit impaired function in the frontal–parietal network (Fassbender & Schweitzer, 2006). The enrichment patterns revealed in this study are consistent with these observations and provide independent evidence at the genetic level emphasizing the roles of the somatomotor network in ADHD.

Our results indicated that SCZ‐associated SNPs were enriched in the frontal–parietal network. Both structural and functional connectivity studies demonstrated that SCZ is associated with severe deficits in the frontal–parietal regions (Pettersson‐Yeo, Allen, Benetti, McGuire, & Mechelli, 2011; van den Heuvel, Mandl, Stam, Kahn, & Hulshoff Pol, 2010). Furthermore, a meta‐analysis of functional neuroimaging studies have reported that during executive function tasks, SCZ patients showed reduced activity in the dorsolateral prefrontal cortex (DLPFC; Minzenberg, Laird, Thelen, Carter, & Glahn, 2009).

One interesting network associated with the set of SCZ‐associated SNPs was the visual regions. SCZ patients show substantial deficits in visual processing (Butler, Silverstein, & Dakin, 2008; Green, Lee, Wynn, & Mathis, 2011). Furthermore, previous studies have reported dysfunctions of the visual cortex in SCZ patients (Cavus et al., 2012; Seymour et al., 2013; van de Ven, Rotarska Jagiela, Oertel‐Knochel, & Linden, 2017). Our results indicated that there is a shared genetic basis between SCZ and the visual functional network. Future studies should investigate more systematically the pathophysiology underlying this relationship.

Like many other studies, the current study used the additive effect model in heritability estimates. Two potential limitations should be noted. First, the additive model examined the narrow‐sense heritability (h 2), which is the proportion of phenotype variance in a population that is attributable to additive genetic variation. It does not consider dominance and epistatic genetic effects; however, these nonadditive genetic effects are believed to be much smaller than the additive genetic effects. A recent meta‐analysis based on 50 years of twin studies, including nearly all published twin studies of complex traits, indicated that 69% of twin studies supported a purely additive genetic model (Polderman et al., 2015). Furthermore, Zhu et al. have applied the extended GCTA approach to estimate dominance genetic variance of 79 quantitative traits using a large sample size of unrelated individuals (~7,000). On average, the dominance genetic variation explained ~3% of the phenotypic variation, which was only a fifth of the additive genetic variation (~15%) (Zhu et al., 2015). Second, the narrow‐sense heritability (h 2), by definition, is non‐negative. Nevertheless, negative values could be observed due to the effect of random errors on very small true heritability. Some studies truncated negative values to zero or did not include them; such a practice might lead to biased estimates of heritability, especially for moderate sample sizes (Wang et al., 1992). To get an unbiased estimate of h 2, the current study allowed for negative values (i.e., did not set heritability estimates between zero and one), even though they were very rare (for a recent in‐depth discussion of negative heritability estimates, see Steinsaltz, Dahl, and Wachter (2018)).

Several lines of future research can further advance our understanding of the genetic architecture of functional connectivity. First, future studies could apply our approach to other big datasets, like the UK biobank, and compare the results with those based on the Chinese samples. Future studies can also combine samples from many different cohorts and conduct comprehensive analyses, such as genome‐wide association analysis, heritability enrichment partitioning analysis, and LD‐MAF stratified heritability analysis (Yang et al., 2015), to examine the contribution of rare genetic variations, which could further unveil the genetic architecture of the multidimensional functional networks. Second, future studies can use longer scan times for the resting‐state scan, which might better characterize individual brain connectivity profiles (D. Wang, Buckner, et al., 2015). Third, the current study used a Chinese sample, whereas the gene sets used in the current study were mainly based on a meta‐analysis of Western participants. Although some studies have reported no heterogeneity between Chinese and European ancestry cohorts (Demontis et al., 2018; Z. Li et al., 2017; Liu et al., 2015), future studies should further replicate our results using SNPs selected based on Chinese populations. Finally, our approaches could be applied to other phenotypes, such as psychiatric disorders and related cognitive functions.

In summary, with a large sample of Han Chinese young adults, the current study revealed for the first time that additive genomic effects explained a considerable proportion of variation in complex brain connectivity profiles. Furthermore, by analyzing the enrichment patterns of gene sets in multidimensional connectivity patterns, the current study provided novel insights regarding the genetic architecture of connectivity profiles in several brain networks, which should help us to better understand the genetic and cognitive mechanisms underlying various psychiatric conditions and eventually guide the development of effective diagnoses and treatments.

CONFLICT OF INTEREST

The authors declare no competing financial interests.

FIGURE S1 Principal component analysis (PCA) of ancestry. The PCA is based on ~3.3 million common SNPs on 1000Genome phase3 data. EAS, East Asian; EUR, European; AFR, African; AMR, Ad Mixed American; SAS, South Asian, StudySample: 1,704 participants used in the current study

TABLE S1. Descriptive summary of the preprocessed genotype data

Click here for additional data file.