The Genetic Architectures of Functional and Structural Connectivity Properties within Cerebral Resting-State Networks

Abstract Functional connectivity within resting-state networks (RSN-FC) is vital for cognitive functioning. RSN-FC is heritable and partially translates to the anatomic architecture of white matter, but the genetic component of structural connections of RSNs (RSN-SC) and their potential genetic overlap with RSN-FC remain unknown. Here, we perform genome-wide association studies (Ndiscovery = 24,336; Nreplication = 3412) and annotation on RSN-SC and RSN-FC. We identify genes for visual network-SC that are involved in axon guidance and synaptic functioning. Genetic variation in RSN-FC impacts biological processes relevant to brain disorders that previously were only phenotypically associated with RSN-FC alterations. Correlations of the genetic components of RSNs are mostly observed within the functional domain, whereas less overlap is observed within the structural domain and between the functional and structural domains. This study advances the understanding of the complex functional organization of the brain and its structural underpinnings from a genetics viewpoint.


Introduction
Structural (SC) and functional connectivity (FC) are vital for healthy cognitive behavior (Buckholtz and Meyer-Lindenberg, 2012; de Lange et al., 2019). Brain regions that show temporally synchronized activity form functionally specialized resting-state networks (RSNs; Yeo et al., 2011), including primary networks (such as the visual or somatomotor network) and higher-order cognitive networks (such as the frontoparietal network, salience network, or default mode network; Bressler and Menon, 2010). Many psychiatric and neurologic disorders have been associated with disruptions within specific RSNs (van den Heuvel and Sporns, 2019). Improving our understanding of the biological principles underlying the concept of structural and functional connectivity within RSNs (RSN-SC/FC) could help elucidate the neural basis of human cognition and disorders associated with disruptions in brain connectivity.
Studies have shown that genetic factors significantly contribute to RSN functional connectivity (twin-based heritability H 2 = 20-40%; Meda et al., 2014;Ge et al., 2017;Adhikari et al., 2018;Miranda-Dominguez et al., 2018;Anderson et al., 2021;Barber et al., 2021). Genomewide association studies (GWAS) on functional connectivity graph measures (Foo et al., 2021) and extrinsic and intrinsic functional organization of RSNs (Zhao et al., 2022) have identified the first genetic variants and sets of genes that make up this genetic component (heritability based on additive common genetic variants, mean h 2 SNP = 13.3%). The genetic determinants of functional connectivity overlap with those of psychiatric disorders (Roelfs et al., 2022). Although RSNs were traditionally discovered based on functional correlation patterns between regions, structural connectivity properties of white matter between the respective brain regions correlate with functional connectivity (Batista-García-Ramó and Fernández-Verdecia, 2018; Mollink et al., 2019;van den Heuvel et al., 2009) to varying degrees across RSNs (Gu et al., 2021). The genetic architecture of RSN structural connectivity has not been investigated to date, but the substantial heritability of multiple properties of major white matter tracts (mean h 2 SNP 25.18-34.9%; Smith et al., 2021;Zhao et al., 2021;Sha et al., 2023) suggests the importance of genetic factors for the anatomic backbone of RSNs. Describing the genetic architecture of both structural and functional RSN connectivity properties as well as annotation and interpretation of the genetic signal can give insight into a biological substrate relevant to a wide variety of neurologic and psychiatric disorders (de Lange et al., 2019) and additionally enables us to estimate to which degree structural and functional RSN connectivity properties are based on a shared genetic source.
This study is aimed to characterize the genetic architecture of within-RSN connectivity properties, both structurally and functionally. 336 and N structural = 23,985; replication N functional = 3408 and N structural = 3412) GWAS are performed on the functional connectivity strength within seven well-known RSNs (Yeo et al., 2011) and a structural connectivity property (fractional anisotropy) within the same RSN definition. These 7 RSN are often used and providing GWAS summary statistics based on the same definition as most neuroimaging studies has the benefit of comparing genetic findings with phenotypic findings. A more granular definition (for example, 17 RSNs by Yeo et al., 2011) was not feasible within our univariate GWAS design, given that the accompanying multiple testing burden would drastically decrease our statistical power to identify and describe genome-wide significant loci (for a multivariate GWAS approach on 17 RSN, see Roelfs et al., 2022). With the polygenic signal from our GWAS, we estimate and partition the heritability and examine the convergence onto genes and biological pathways, with the purpose of aiding the biological interpretation of results and providing meaningful starting points for functional follow-up experiments (Uffelmann and Posthuma, 2020). We examine genetic correlations between different RSNs, as well as across structural and functional domains. These genetic correlation analyses are extended to the locus level to facilitate the prioritization of possible pleiotropic loci for future studies (Werme et al., 2022). Altogether, we focus on the translation of RSN-associated genetic loci into biological interpretation and provide insights into the genetic architecture of within-RSN functional and structural connectivity properties.

Materials and Methods
A flowchart that describes all methods used in this manuscript is displayed in Extended Data Figure 1-3.

Sample
The UK Biobank (UKB) is a resource with genomic and imaging data of volunteer participants (Sudlow et al., 2015). The National Research Ethics Service Committee North West-Haydock ethically approved this initiative (reference 11/NW/0382) and data were accessed under application #16406. Combined single nucleotide polymorphism (SNP)-genotypes and neuroimaging data of N = 40,682 participants have been available since January 2020. From all new subjects in the latest neuroimaging release (January 2020), we randomly assigned 5000 subjects to a holdout set for validation. Subsetting the total sample to subjects with all neuroimaging data necessary to construct our phenotypes as described below, resulted in N functional = 37,017 and N structural = 36,645. We only included subjects for which the projected ancestry principal component score was closest to and ,6 SD from the average principal component score of the European 1000 Genomes sample based on Mahalanobis distance. This procedure has been described in previous publications by our group (Jansen et al., 2020) and the number of non-European exclusions are displayed in Extended Data Figure 1-2. Other exclusion criteria were withdrawn consent, UKB-provided relatedness, discordant sex, or sex aneuploidy (Extended Data Fig. 1-2). Further quality control on genomic and neuroimaging data are described below and resulted in the sample sizes and sample characteristics as displayed in Extended Data Figure 1-1.

Genotype data
The genotype data used in this study for discovery and validation were obtained from the UK Biobank Axiom and the UK BiLEVE Axiom arrays. These Affymetrix arrays cover 812,428 unique genetic markers and overlap 95% in SNP content. This number of SNPs was increased to 92,693,895 by imputation conducted by UKB. Variants were imputed using the Haplotype Reference Consortium and the UK10K haplotype panel as reference. We applied our in-house quality control pipeline in addition to quality control performed by UKB. This procedure excluded SNPs with low imputation scores (INFO , 0.9), low minor allele frequency (MAF , 0.005) or high missingness (.0.05), multiallelic SNPs, indels, and SNPs without unique rs-identifiers. A total of 9,380,668 SNPs passed quality control and were converted to hard call SNPs using a certainty threshold of 0.9 for further analyses.

Neuroimaging data Preprocessing and connectome reconstruction
The UKB scanning protocol and processing pipeline is described in the UKB Brain Imaging Documentation . For this study, we made use of the available resting-state functional brain images (rsfMRI) and multiband diffusion brain images (DWI) together with T1 surface model files and structural segmentation from FreeSurfer (Fischl, 2012). These three types of data were used as input for the structural and functional pipeline of CATO (Connectivity Analysis TOolbox; de Lange et al., 2022). Before this, UKB performed preprocessing on DWI and rsfMRI data as described in the UKB Brain Imaging Documentation .
In CATO's structural pipeline, additional preprocessing of DWI files was performed in FSL (Jenkinson, et al., 2012) by computing a DWI reference image based on the corrected diffusion-unweighted (b0) volumes, computing the registration matrix between DWI reference image and the anatomic T1 image, and registering the Freesurfer segmentation to the DWI reference image. The surface was parcellated based on the Cammoun subparcellations of the Desikan-Killiany atlas including 250 cortical regions (Cammoun et al., 2012). We reconstructed the diffusion signal with diffusion tensor imaging (DTI), a deterministic method that is robust and relatively simple compared with more advanced diffusion reconstruction methods (de Lange et al., 2022). In CATO, the fiber assignment by continuous tracking (FACT) algorithm (Mori et al., 1999) is used to reconstruct fibers and fractional anisotropy (FA) was used as weights of reconstructed fibers. FA is a relatively simple measure of white matter integrity and correlates with axon density, size and myelination (Beaulieu, 2002). The structural connectivity matrix was built out of all fiber segments that connected two regions in the atlas. Additional filters were applied, namely a minimal FA of 0.1, minimal length of 30 mm and having 2 or more number of streamlines.
The functional pipeline in CATO consisted of similar steps. First, we computed an rsfMRI reference image by averaging all rsfMRI frames in FSL and subsequently registered this reference image and the T1 image in FreeSurfer. Second, we parcellated the surface based on the same atlas as in the structural pipeline (to enable structure-function comparison in downstream analyses) and we registered the T1 parcellation to the rsfMRI image. Third, motion metrics were estimated, and time-series were corrected for covariates (linear trends and first order drifts of motion parameters and the mean signal intensity of voxels in white matter and cerebrospinal fluid and of all voxels in the brain) by regression. Fourth, time-series were passed through bandpass filtering (frequencies 0.01-0.1) and scrubbing (max FD = 0.25, max DVARS = 1.5, min violations = 2, backward neighbors = 1, forward neighbors = 0). Fifth, the functional connectivity matrix was computed by the Pearson's correlation coefficient of the average signal intensity of every pair of brain regions across the frames that survived filtering.

Quality control
The UKB scanning and preprocessing protocol includes filters for outliers based on manual QC and an advanced classifier described elsewhere (Alfaro-Almagro et al., 2018). We excluded a small number of subjects that UKB identified as outliers and placed in an "unusable" folder. The UKB main documentation  suggests a second set of UKB data fields that can be used as outlier criteria. Outlier subjects are defined as subjects that score for any of the values .3 interquartile ranges above the upper quartile or below the lower quartile. Outlier criteria included measures that describe the discrepancy between the T1-weighted, rsfMRI and DWI images and the population average template after LINEAR and NON-LINEAR alignment, the amount of nonlinear warping necessary to map a subject to the standard template, the signal-to-noise ratio in rsfMRI, the mean rfMRI head motion averaged across space and time points and the total number of outlier slices in DWI volumes. We extended this recommended list with connectome specific measures, including the average prevalence of all connections present and absent in the reconstructed brain network of a subject (low average prevalence scores indicate the presence of odd connections and high values indicate the absence of common connections), the sum of numberof-streamlines and average FA of all connections in the reconstructed brain network of a subject. The number of exclusions can be viewed in Extended Data Figure 1-2.

Phenotype reconstruction
In this study, the phenotypes of interest were functional and structural connectivity properties (FC; SC) within seven resting-state networks (RSNs) that previously have been identified (Yeo et al., 2011) and are commonly used in (clinical) neuroimaging studies: the default mode network, ventral attention network, dorsal attention network, visual network, limbic network, somatomotor network and frontoparietal network. Each of the 250 cortical regions of the reconstructed structural and functional connectomes were assigned the ratio to what extent they belonged to each of these seven RSNs, using a mask created and validated elsewhere (see Wei et al., 2019). Each connection was then weighted by multiplying the ratios of the two regions involved in the particular RSN. FC and SC properties within the RSNs were, respectively, calculated as the mean correlation and mean fractional anisotropy of the connections within the RSN. We also computed two global FC and SC phenotypes as the mean correlation and mean fractional anisotropy of all available connections, to be able to correct for connectivity that is nonspecific to RSNs in downstream analyses.

Statistical analyses SNP-based GWAS
To identify common genetic variants involved in FC within each of the seven RSN, we performed seven single nucleotide polymorphism (SNP)-based genome-wide association studies (GWAS) in PLINK2 (Purcell et al., 2007). Also, for the SC property within each of the seven RSN, a SNP-based GWAS was performed. It is common practice to include a global FC or SC estimate as covariate in GWAS to capture associations that are driven by the level of connectivity within an RSN regardless of the level of connectivity throughout the whole brain. It has become apparent that this risks the introduction of collider bias (inducing false-positives; Day et al., 2016). Here, we build on recent developments in statistical genetics that have provided multiple methods that allow for post-GWAS analyses conditional on global connectivity. Therefore, we used the global FC and global SC phenotypes to run two additional SNP-based GWAS, for which the summary statistics were used for conditional downstream analyses. The total number of GWASs was therefore 16. In order to correct for population stratification during GWAS, a principal component analysis was performed in FlashPCA2 (Abraham et al., 2017) using only independent (r 2 , 0.1), common (MAF . 0.01), and genotyped SNPs or SNPs with very high imputation quality (INFO = 1). The first 30 principal components were used as covariates in all GWASs, together with sex, age, genotype array, Townsend deprivation index (a proxy of socio-economic status), and general neuroimaging confounders as well as FC/SC specific covariates (recommended by Alfaro-Almagro et al., 2020). The general set included handedness, scanning site, the use of T2 FLAIR in Freesurfer processing, intensity scaling of T1, intensity scaling of T2 FLAIR, scanner lateral (X), transverse (Y), and longitudinal (Z) brain position, and z-coordinate of the coil within the scanner. FC-specific and SC-specific covariates were, respectively, intensity scaling and echo time of rsfMRI, and intensity scaling of DWI. For reasons of collinearity, we ran principal component analysis on all covariates (excluding the population stratification principal components) and retained those principal components that explained . 99% of variance. Rare variants (MAF , 0.005) and SNPs with high missingness (.5%) were excluded from GWAS and male X variants were counted as 0/1. The genome-wide significance threshold was a = (0.05/1,000,000/16 =) 3.13 Â 10 À9 according to the Bonferroni correction for multiple testing.

SNP-based heritability
SNP-based (h 2 SNP ; or narrow-sense) heritability represents the proportion of phenotypic variance that can be explained by common additive variation. In contrast, broadsense heritability captures the total genetic contribution to the phenotype and is often based on family studies (Visscher et al., 2008). We applied linkage disequilibrium score regression (LDSC) on the SNP-based GWAS summary statistics of all 16 phenotypes to estimate h 2 SNP using precomputed LD scores from 1000 Genomes EUR, as provided by the LDSC developers.
For the phenotypes with enough polygenic signal to run LDSC (l . 1.02), we investigated whether certain functional categories in the human genome were enriched for h 2 SNP . The ratio of the proportion of h 2 SNP in a certain category to the proportion of SNPs in the category equaled the enrichment value. We corrected the level of significance for multiple testing to a = (0.05/28)/11 = 1.62 Â 10 À4 .

Functional annotation
Functional mapping and annotation (FUMA) is a webbased platform that can be used to functionally map and annotate SNPs that appear significant during GWAS. We uploaded summary statistics to FUMA if GWAS identified at least one genome-wide significant SNP. Candidate SNPs were defined as all SNPs in linkage disequilibrium (LD) r 2 . 0.6 with an independent genome-wide significant SNPs (r 2 , 0.6). Annotation was subsequently performed using ANNOVAR (K. Wang et al., 2010), RegulomeDB score (Boyle et al., 2012) and ChromHMM (Ernst and Kellis, 2012). Lead SNPs were defined as independent SNPs r 2 , 0.1. Genomic loci were constructed by taking all independent significant SNPs r 2 , 0.1 with LD blocks within 250 kb distance and independent significant SNPs r 2 ! 0.1. Within every locus, SNPs were mapped to genes using three methods: positional mapping, expression quantitative trait loci (eQTL) mapping or chromatin interaction mapping. SNPs were positionally mapped to genes if their physical distance was ,10 kb. Mapping based on eQTLs relied on known associations between SNPs and the gene expression of genes within a 1-Mb window, from BRAINEAC frontal, occipital, temporal and cerebral cortex (Ramasamy et al., 2014), GTEx v8 cerebral cortex (GTEx Consortium, 2020) and xQTLServer dorsolateral prefrontal cortex (Ng et al., 2017). Chromatin interaction mapping was based on established 3D DNA-DNA interactions between SNP and gene regions from Hi-C databases in cortex tissue [PsychENCODE by D. Wang et al. (2018), Giusti-Rodriguez et al. (2019) and GSE87112 (Schmitt et al., 2016)]. To restrict chromatin interaction mapping to plausible biological interactions, we only included interactions where one region overlapped with an enhancer [as predicted by the Roadmap Epigenomics Consortium (2015) in cortex tissue] and the other region overlapped with a promoter (250 bp upstream to 500 bp downstream of the transcription start site as well as predicted by the Roadmap Epigenomics project in cortex tissue). A false discovery rate (FDR) threshold of 1 Â 10 À5 was used, as recommended in previous literature (Schmitt et al., 2016).

Statistical fine-mapping of genome-wide significant loci
Statistical fine-mapping methods can aid in determining the probability of variants identified through GWAS being causal. We applied FINEMAP to the genome-wide significant loci defined by FUMA, which is a Bayesian statistical fine-mapping tool that estimates the posterior probability of a specific model by combining the prior probability and the likelihood of the observed summary statistics (Benner et al., 2016). The posterior probabilities are used to calculate the posterior inclusion probability (PIP) of a SNP in a model and the minimum set of SNPs needed to capture the SNPs that most likely cause the association (Schaid et al., 2018). We set the maximum number of causal SNPs to k = 10 and report on the 95% credible set of the model with k causal SNPs reaching the highest probability. LDstore (Benner et al., 2017) was used to estimate the pairwise LD matrix of SNPs from quality controlled genomic data of the discovery sample. Only SNPs with PIP . 0.10 were used for interpretation in the main text, but all are reported on in Extended Data Figure 3-2.

Gene-based GWAS
Performing GWAS on the level of genes has been suggested to be more powerful than GWAS on the level of SNPs (de Leeuw et al., 2015). Therefore, the 16 SNPbased GWAS summary statistics were used to perform 16 gene-based GWAS in MAGMA (multimarker analysis of genomic annotation) v1.08 (de Leeuw et al., 2015). A mean SNP-wise model was applied (with the UKB European population serving as an ancestry reference group) to test the joint association of all SNPs within 18,850 genes with FC/SC within RSNs. The genome-wide significance threshold was adjusted for multiple testing to a = (0.05/18,850)/ 16 = 1.66 Â 10 À7 .

Gene-set analysis
We set out to prioritize the associations from genebased GWAS and gain more insight into the biological pathways associated to RSN-FC/RSN-SC. In order to identify gene-sets specific for RSN-FC/RSN-SC, we ran conditional gene-set analyses in MAGMA conditioning on the global FC or SC, respectively. Pathways were represented as gene-sets from Gene Ontology (GO) molecular functions, cellular components and biological processes, and curated gene-sets from MsigDB v7.0 (sets C2 and C5; Liberzon et al., 2011). Protein-coding genes served as background genes. The threshold for gene-sets reaching significance was corrected for multiple testing to a = (0.05/7246)/14 = 4.93 Â 10 À7 .

Genome-wide genetic correlations
We aimed to test whether genetic signals were generally similar across the FC/SC within various RSNs. To assess the overlap in genetic architecture between FC/SC within RSNs while taking the influence of global FC/SC into account, we designed a genetic correlation (r g ) analysis pipeline. This pipeline consisted of three steps. (1) In the first step, genome-wide r g between 42 combinations of RSNs were estimated using LDSC [a = (0.05/42=) 1.19 Â 10 À3 ]. The summary statistics of SNP-based GWAS were used as input for LDSC. We excluded FC within the visual network, because both the l (,1.02) and ratio (.0.20) values were out of bound for LDSC.
(2) For all RSNs included in a significant bivariate r g , additional r g with global FC/SC were calculated in LDSC.
(3) If one or both RSNs from the significant bivariate r g showed additional significant r g with global FC/SC, we recalculated of the genome-wide r g between the two RSNs with global FC/SC taken into account. Since such residual genome-wide r g analyses are not implemented in LDSC, we applied genomic structural equation modelling (genomic SEM; Grotzinger et al., 2019). Genomic SEM is a method that enables to model the multivariate genetic architecture and covariance structure of complex traits using GWAS summary statistics and allows for sample overlap. We modelled residual covariance between RSNs as the covariance between the residual variance of the two RSNs involved after taking the global factor into account (Fig. 1). The model was then fit using diagonally weighted least square estimation.

Local genetic correlations
The genome-wide r g' s described above are an average correlation of genetic effects across the genome, implicating that contrasting local r g' s are possibly cancelling each other out. Running r g analysis on a locus level has the potential to uncover loci that show genetic similarity between traits. For this purpose, we adopted a three-step local r g analysis pipeline similar to the genome-wide r g analysis approach described above. All three steps were performed in LAVA (Werme et al., 2022), a local r g analysis R package, using SNP-based GWAS summary statistics as input. We followed the suggested sample overlap procedure (as described on https://github.com/josefin-werme/LAVA) to enable LAVA to model shared variance because of sample overlap as residual covariance and consequently remove upward bias in local r g estimates (Werme et al., 2022). Since our GWASs included European samples, the 1000 Genomes Figure 1. Path diagram of genomic SEM model. The summary statistics of two RSNs that have shown to significantly correlate with global connectivity will be used as input together with summary statistics of the global connectivity GWAS. In this way, r g between the unique components (u) of the two RSNs can be estimated while taking global connectivity into account. The example in this diagram shows the global and unique genetic effects on functional connectivity (FC) for the ventral attention network (VAN) and default mode network (DMN), but a similar model was used for other RSN pairs and for measures of structural connectivity (SC). This method is embedded in a flowchart that describes the sample (see Extended Data Figs. 1-1 and 1-2 for sample characteristics and exclusion criteria) and all methods used in this manuscript (Extended Data Fig. 1-3).
Phase 3 European data served as genotype reference and formed the basis of the locus definition file. For every locus, the first step of our pipeline consisted of estimating local bivariate r g between 49 combinations of RSNs. However, RSNs that were devoid of heritable signal (p . 1 Â 10 À4 ) in the locus were excluded from local bivariate r g analysis to ensure interpretability and reliability. A total of 774 bivariate tests were performed across 337 loci, leading to an adjusted significance threshold of a = (0.05/774=) 6.46 Â 10 À5 . In the second step, RSNs that showed significant local r g were additionally tested for r g in that locus with global FC/SC. Note that if this was not possible, because global FC/SC showed no significant heritability in that locus, the local bivariate r g between RSNs could not be biased by global FC/SC. If one or both RSNs did show additional significant r g with global FC/SC, we ran a partial local r g between the RSNs conditioned on the SC-global and/or FC-global phenotype in step three. If the partial local r g between the RSNs no longer remained significant, we concluded that the initial r g was driven by global FC/SC and did not reflect genetic overlap specific for these RSNs.

Polygenic score
The variance explained in RSN-FC/RSN-SC by our GWAS findings was investigated to test the robustness of our findings using polygenic score (PGS) estimation in PRSice-2 (Choi et al., 2019). We applied a two-phase approach to obtain p-values unaffected by overfitting and therefore split our holdout sample in a target (N = 1818) and validation set (N = 1824). In phase 1, SNP-based summary statistics (MAF . 0.1, chromosome X excluded) from the discovery sample together with genotype data of the target set were used to find the optimal p-value threshold. PRSice-2 uses high-resolution thresholding and clumping of genotype data, and we included the same covariates as during discovery GWAS. In phase 2, the model with the best fit from phase 1 was applied on the validation set to obtain the variance explained (R 2 ).

Replication of lead SNPs
In the design of this study, a hold-out sample of N FC = 3408 and N SC = 3412 was reserved for PGS analysis. We applied an earlier described method (see Okbay et al., 2016) to internally validate our discovery lead SNPs. This formula describes the probability of a discovery lead SNP i being significant in a replication sample as with a representing an a level of 0.05, U the cumulative normal distribution function, U À1 the inverse normal distribution function, s rep;i the standard error of SNP i in the replication GWAS and b i the winner's curse adjusted association estimate of SNP i. Winner's curse is the occurrence of overestimated effect sizes that are induced by significance thresholding (Palmer and Pe'er, 2017). We applied winner's curse correction using the mean of the normalized conditional likelihood (Ghosh et al., 2008)

Data availability
Genome-wide summary statistics will be made publicly available via https://ctg.cncr.nl/software/summary_statistics/ upon publication.

GWAS of RSN-SC and RSN-FC properties identify six genome-wide significant loci
Following previously described procedures , we started our analysis by grouping cortical areas into seven RSN as defined by Yeo et al. (2011;visual, somatomotor, limbic, dorsal attention, ventral attention, frontoparietal, and default-mode network; Extended Data Fig. 1-3) and calculating the mean functional and structural connectivity within the RSNs in UK Biobank subjects (discovery N FC = 24,336 and N SC = 23,985; replication N FC = 3408 and N SC = 3412). Within-RSN functional connectivity strength (from now on referred to as RSN-FC) was measured as the average correlation between the activation signals of brain regions within each RSN over time. A property of within-RSN structural connections (from now on RSN-SC) was measured as the average fractional anisotropy (FA) of white matter connections between brain regions within each RSN (see Materials and Methods). FA values are believed to reflect a metric of efficiency of or capacity for information transfer across white matter pathways and are sensitive to myelin content (Beaulieu, 2002). Discovery GWAS were performed for the FC and SC within every RSN and identified 518 genome-wide significant variants (or single nucleotide polymorphisms; SNPs) at p , (5 Â 10 À8 /16 =) 3.13 Â 10 À9 located in six genomic loci (Fig. 2): three for visual network-SC, one for limbic network-FC, and a shared locus for frontoparietal network-FC and somatomotor network-FC (Extended Data Fig. 2-6). These loci do not seem simply driven by overall connectivity properties, given that none of these six loci showed a genome-wide significant association with global FC or SC.
The proportions of phenotypic variance explained by additive genetic effects of GWAS SNPs (SNP-based heritability; h 2 SNP ) were moderately higher for RSN-SC (M = 13.59%, SD = 1.79%) than those observed for RSN-FC (M = 6.71%, SD = 3.36%; Extended Data Fig.  2-8). We did not find evidence for enrichments of h 2 SNP in functional genomic categories after Bonferroni-correction (Extended Data Fig. 2-9). The linkage disequilibrium score regression (LDSC) intercept approached one for all phenotypes, indicating limited bias from population stratification.

Axon guidance and synaptic functioning genes implicated in visual network-SC GWAS
We continued by examining the possible functional consequences of the SNPs involved in RSN-FC and RSN-SC. SNPs in linkage disequilibrium (i.e., in correlation; LD; r 2 ! 0.6) with the Bonferroni-corrected genome-wide significant SNPs from the GWAS which also had suggestive p-values (,1 Â 10 À5 ) and a minor allele frequency (MAF) . 0.005 were defined as candidate SNPs (Extended Data Fig. 2-7). We subsequently mapped candidate SNPs to genes using three strategies (Materials and Methods): positional mapping if the SNP was within ,10 kb of a gene, expression quantitative trait loci (eQTL) mapping if the SNP was known to affect expression of the gene within 1 Mb in cortex, and chromatin interaction mapping if there was a significant Hi-C interaction between the SNP and a nearby or distant gene in cortex tissue (Extended Data Fig. 3-1).
For visual network-SC, an exonic nonsynonymous SNP located in exon 1 of AC007382.1 (rs711244, p = 1.42 Â 10 À12 ) was among the candidate SNPs in the locus on chromosome 2. The function of AC007382.1 is unknown, but it has been associated with amygdala volume previously (Mufford et al., 2021). A SNP 5 kb from AC007382.1 (rs35050623) was also among the most likely causal SNPs as identified by fine-mapping the locus (posterior inclusion probability; PIP = 0.12). The loci on chromosome 10 and 7 were "unsolved" by FINEMAP (all PIPs , 0.1), which can occur because of a combination of extensive linkage disequilibrium and small effect sizes (Extended Data Fig. 3-2). Exonic synonymous SNPs were found in exon 7 of FAM175B and exon 12 of SEMA3A in the loci on chromosome 10 and 7, respectively. The transcript of FAM175B is a component of an enzyme complex that deubiquitinates Lys-63 linked chains to control protein function (Cooper et al., 2010). Experimental studies have suggested that such deubiquitination can regulate synaptic transmission and synaptic plasticity (DiAntonio and Hicke, 2004). SEMA3A contained multiple intronic SNPs associated with visual network-SC with high combined annotation dependent depletion (CADD) scores (11 SNPs with CADD . 12.37), which are usually considered reducing organismal fitness and correlating with molecular functionality and pathogenicity (Kircher et al., 2014). The product of SEMA3A is known as a key regulator of axon outgrowth during the establishment of correct pathways in the developing nervous system (Zhou et al., 2008).
We additionally mapped 46 visual network-SC candidate SNPs to METTL10, because of their established associations with METTL10 expression levels in fetal and adult cerebral cortex tissue (eQTL mapping) as well as their chromatin interaction. METTL10 encodes a methyltransferase that catalyses the trimethylation of eukaryotic elongation factor a (eEF1A) at Lys-318, a key regulator of ribosomal translation (Jakobsson et al., 2018). Visual network-SC SNPs were also mapped to the METTL10-FAM53B readthrough (RP11-12J10.3) and FAM53B gene, because of known chromatin interaction in fetal and adult cerebral cortex tissue (Fig. 3a). FAM53B is required for Wnt signaling, a pathway important for cell regeneration (Kizil et al., 2014). Lastly, positional mapping of candidate SNPs within a 10-kb window of a gene resulted in the identification of VIT, STRN, and HEATR5B genes for visual network-SC (Extended Data Fig. 3-1). Two intronic SNPs within STRN (rs2003585, rs2691112) were also included in the 95% credible causal SNP set during finemapping with a posterior probability of being causal (PIP) of 0.14 and 0.12, respectively (see Extended Data Fig. 3-2 for all SNPS in the 95% credible sets identified by finemapping).

Annotation of identified loci for RSN-FC
We observed two exonic nonsynonymous SNPs in exon 19 (rs2274224, p = 1.771 Â 10 À10 ) and exon 25 (rs2274223, p = 1.22 Â 10 À5 ) of the PLCE1 gene to be associated with limbic network-FC. The PLCE1 gene encodes for the phospholipase C e 1, which mediates the production of two second messengers that regulate cell growth, differentiation, and gene expression (Rao et al., 2017). The high CADD scores (17.35 and 17.48, respectively) suggest deleteriousness of these two exonic nonsynonymous SNPs. rs2274224 was also among the FINEMAP 95% credible set of six SNPs Research Article: New Research within PLCE1, although the intronic rs10786156 had the highest probability of being causal (PIP = 0.36). Additionally, four intergenic SNPs within the same locus were located near the NOC3L gene.
On chromosome 10, eight SNPs associated with limbic network-FC were mapped to the CYP2C8 gene based on eQTL mapping (Fig. 3b). Expression of CYP2C8 results in an enzyme important for drug metabolism (Backman et al., 2016). One of CYP2C8 substrates, the nonselective monoamine oxidase inhibitor phenelzine, is known to target the nervous system and is clinically prescribed as treatment for major depressive disorder (Q. Wang et al., 2019). A large body of research has verified the association between major depressive disorder and changes in limbic network functional connectivity, as well as with other RSNs (for a meta-analysis, see Kaiser et al., 2015).
The annotation of SNPs in the locus that was shared between frontoparietal and somatomotor network-FC revealed only intergenic candidate SNPs (enrichment = 2.15, p = 5.09 Â 10 À9 ), which convolutes biological interpretation but is a common observation for complex traits . The nearest genes to the candidate SNPs in this locus were PAX8 and IGKV1OR2-108 (respectively, 29 and 53 kilobase distance). The likely deleterious SNP rs199993536 (CADD = 19.87) closest to PAX8 was also the most likely causal SNP in the 95% credible set for both frontoparietal network-FC (PIP = 0.19) and somatomotor network-FC (PIP = 0.42), with a probability of 0.08 that this is the shared causal variant between the two networks. PAX8 encodes a transcription factor that is considered to regulate the expression of genes important for thyroid development (Blake and Ziman, 2014) and the production of thyroid hormone (Di Magliano et al., 2000). FC within both the somatomotor and frontoparietal network is reduced in individuals with subclinical (Kumar et al., 2018) and clinical hypothyroidism (Singh et al., 2015).

Default mode network-FC genes associated with Alzheimer's disease
We next performed gene-based GWAS for the FC and SC within every RSN (Fig. 4). We detected two Bonferronicorrected genome-wide significant genes additional to the mapped genes described above by combining association statistics from neighboring variants within a single gene (Extended Data Fig. 4-1). Visual network-FC was associated with APOC1 (z = 5.15, p = 1.31 Â 10 À7 ), and for default mode network-FC APOE was found to be associated (z = 5.13, p = 1.43 Â 10 À7 ). APOC1 and APOE are both located within the 19q13.2 locus and are well-known risk factors for Alzheimer's disease (Emrani et al., 2020).
In order to determine whether there is genetic overlap between Alzheimer's disease (Wightman et al., 2021) and default mode network-FC, we performed local genetic correlation (r g ) analysis (see Materials and Methods; Extended Data Fig. 4-2). For default mode network-FC, we detected two loci on chromosome 12 (BP 64,403,114,643) and 19 (BP 45,040,893,307) which showed significant local r g at p , (0.05/71=) 7.04 Â 10 À4 with Alzheimer's disease (Extended Data Fig. 4-3). Given the negligible heritability of global FC in these loci (univariate p = 0.27 and p = 0.01, respectively, whereas p = 1.30 Â 10 À5 and p = 1.62 Â 10 À8 for default mode network-FC) we conclude that these local genetic associations with Alzheimer's disease are not driven by total brain connectivity. The locus on chromosome 12 showed a positive r g (r ) between Alzheimer's disease and default mode network-FC (BP 64,403,114,643, r = 0.69, 95% CI = 0.35-1.00, p = 3.25 Â 10 À4 ). Interestingly, this locus has been identified in a previous GWAS for hippocampal atrophy, a biological marker of Alzheimer's disease (Bis et al., 2012). Negative r g between Alzheimer's disease and default mode network-FC was observed in the locus on chromosome 19 (BP 45,040,893,307, r = À0.56, 95% CI = À0.82 to À0.38, p = 9.23 Â 10 À9 ), indicating that lower default mode network-FC was associated with higher genetic risk of Alzheimer's disease. Note that this larger defined locus showed weak heritability (p = 0.014) for visual network-FC despite the significance of APOC1 in the gene-based GWAS, which would make genetic correlation estimates with Alzheimer's disease unreliable and uninterpretable (Werme et al., 2022). Therefore, Alzheimer's disease seems to show genetic overlap specifically with default mode network-FC in this locus.

Looking for biological pathways through gene-set analysis
We looked for convergence of the genetic signal for RSN-FC/RSN-SC onto 7252 MSigDB (Liberzon et al., 2011) pathways using MAGMA gene-set analysis, a useful method for further functional interpretation (Uffelmann and Posthuma, 2020). We conditioned our analyses on the gene-based GWAS summary statistics for global FC/SC in an effort to capture RSN-specific pathways. Five pathways showed an association with four RSNs after Bonferroni correction for the number of pathways tested per trait (Extended Data Fig. 4-4 displays the associations of all pathways tested for all traits). These included blood vessel morphogenesis (GO, p = 3.30 Â 10 À6 ) and vasculature development (GO, p = 4.94 Â 10 À6 ) pathways for limbic network-SC, the Parkinson's disease pathway (KEGG, p = 2.10 Â 10 À6 ) for somatomotor network-SC, the pathway for positive regulation of mesenchymal cell proliferation (GO, p = 2.64 Â 10 À6 ) in default mode network-FC, and the pathway for negative regulation of histone methylation (GO, p = 5.23 Â 10 À6 ) for the dorsal attention network-FC. These five pathways did however not survive further Bonferroni correction for the number of traits tested [a = (0.05/7246)/ 14 = 4.93 Â 10 À7 ]. Therefore, it cannot be concluded that these biological processes are involved in the genetics of FC and SC within RSNs.

Examining overlap between structure and function per RSN through genetic correlations
As SC strength has been noted to correlate with FC strength on the phenotypic level (Mollink et al., 2019), we sought to investigate the correlations between FC and SC properties within each RSN on a genetic level. Genomewide genetic correlations (r g ) were estimated using SNPbased summary statistics (Fig. 5). We observed no nominally significant genome-wide r g s between SC and FC in any of the RSNs (Extended Data Fig. 5-1). Genomewide r g estimates ranged from À0.19 (SE = 0.15, p = 0.19) in the dorsal attention network and 0.23 (SE = 0.23, p = 0.30) in the frontoparietal network.
Strongly localized or opposing local r g' s possibly may go undetected, since genome-wide r g' s are an average of the shared genetic association signal across the genome. We examined whether such relationships between SC and FC within any given RSN exist by performing local r g analysis (Werme et al., 2022), although we did not identify any significant r g on a locus level either (Extended Data Table 1-1).

Genome-wide and local genetic correlations within the functional and structural domain
We examined the shared genetic signal across RSNs within the same domain by conducting genome-wide r g Figure 5. Global r g (6SE) between functional connectivity (FC) and structural connectivity (SC) within the same RSN. Genetic correlations as performed in LDSC do not show estimates significantly different from 0 (Extended Data Fig. 5-1). Additional estimation of local r g did not yield significant overlapping loci between SC and FC within each RSN either. The colors correspond to the RSN colors in Figure 2 and 4. analyses ( Fig. 6; Extended Data Fig. 5-1). For functional connectivity strength within the default mode network and within the ventral attention network, a significantly shared genetic signal was observed after Bonferroni correction (r g = 0.52, SE = 0.16, p = 1.00 Â 10 À3 ). This association was not driven by global FC as neither default mode network-FC nor ventral attention network-FC were genetically correlated with global FC (r g = 0.19, SE = 0.18, p = 0.29; r g = 0.26, SE = 0.19, p = 0.18, respectively). Note that this positive r g does not imply simultaneous functional activation of these two RSNs or their involvement in similar cognitive tasks (which would contradict previous research, Menon and Uddin, 2010), but suggests that variants that influence default mode network-FC generally tend to influence ventral attention network-FC in the same direction.
For the structural connectivity property, we observed multiple significant genome-wide r g' s (p , 1.19 Â 10 À3 ) across RSNs, although many of these were also correlated with global SC (Extended Data Fig. 5-1). To determine whether genetic overlap of structural connectivity across RSNs could be accounted for by global SC, we used genomic SEM to compute residual r g estimates across RSN-SC while taking global SC into account (see Materials and Methods). As none of the residual r g estimates remained significant, we conclude that global SC likely accounts for the observed genetic overlap across RSN-SC.
We extended our investigation into shared genetic signal across RSNs beyond the global to the local scale. Eighteen loci showed Bonferroni corrected significant r g s when comparing RSNs within the functional domain (Table 1). These were all highly positive (mean r = 0.84, SD = 0.09) and were not confounded by global FC. When comparing RSNs within the structural domain, local r g analysis revealed only one positively correlated locus between dorsal attention network-SC and frontoparietal network-SC (chr15:39238841-40604780, local r g (r ) = 0.85, p = 9.51 Â 10 À7 ; Table 1). A complete overview of local r g results can be found in Extended Data Table 1-1.

Lead SNP validation and polygenic score prediction
We examined the replicability of the discovery lead SNPs as defined in FUMA in our holdout sample (N FC = 3408; N SC = 3412). From these six lead SNPs, we estimated to replicate three (exact 3.99) at a = (0.05/n lead SNPs per trait) in our holdout sample given their winner's curse corrected effect size and the sample sizes of the discovery and replication samples. Observations showed three discovery lead SNPs to be significant (Extended Data Fig. 2-4). Extended Data Figure 2-5 shows the probability distributions for all six discovery lead SNPs to be significant at increasing replication sample sizes.
Additionally, the variance that could be explained in RSN-FC/RSN-SC by polygenic scores based our GWAS associations was considered. In PRSice-2, the summary statistics of the discovery GWAS were used to find the best p-value threshold for polygenic scores in the target set. The application of this optimal prediction model in our validation set explained on average 0.28% and 0.35% of the variance across RSN-FC and RSN-SC, respectively (Extended Data Fig. 2-3). Note that this R 2 value is on average 2.75% (SC) to 6.89% (FC) of the h 2 SNP , which is a comparable with other studies with relatively small sample sizes and is expected to climb close to h 2 SNP with increasing sample sizes (Choi et al., 2020). . Genome-wide r g across RSN measures of (a) functional connectivity (FC) and (b) structural connectivity (SC). If one of the two RSNs showing significant LDSC r g showed additional significant r g with global FC/SC, we instead report the residual r g (r g between the two RSNs while taking global FC/SC into account in Genomic SEM; see Materials and Methods and Fig. 1). The significant r g that survived correction for multiple testing (p , 1.19 Â 10 À3 ) is indicated with an asterisk (*).

Discussion
Mapping the genetic components of resting state networks (RSNs) may provide insight into the etiology of brain function and brain disorders. RSNs are typically defined using functional connectivity (FC) patterns across brain regions, and structural connectivity properties (SC) between these regions correlate with FC across RSNs (Gu et al., 2021). The aim of this study was to gain more insight into the genetic underpinnings of structural and functional connectivity properties (SC; FC) within a framework that respects the brain's hierarchical functional architecture. With the use of GWAS and annotation we observe that genetic variation in RSN-FC (e.g., limbic network-FC and default mode network-FC) impacts biological processes related to brain disorders (major depressive disorder and Alzheimer's disease, respectively) that have previously been associated with FC alterations in those RSN. We further identify genes for visual network-SC that are involved in axon guidance and synaptic functioning. The genetic component of RSNs overlaps mostly within the functional domain, whereas less overlap is observed within the structural domain and between the functional and structural domains.
For FC within RSNs (RSN-FC), we detect biologically interpretable results for the default mode and limbic network-FC. For default mode network-FC, we observe APOE as a genome-wide significant gene. The default mode network is hypothesized to relate to Alzheimer's disease through the role of default mode network-FC in memory consolidation (Fox and Raichle, 2007) and through cortical atrophy spreading over default mode network regions over time (Seeley et al., 2009). Here, we complement earlier phenotypic observations that link Alzheimer's disease to default mode network-FC (Chiesa et al., 2017) by now also showing genetic correlations in two loci between Alzheimer's disease and default mode network-FC. Functional follow-up would be necessary to investigate how the variants and genes in these loci affect default mode network-FC. The limbic network is commonly known for its involvement in emotion regulation, episodic memory, and action-outcome learning (Rolls, 2019) and has been associated with mood disorders such as major depression disorder and bipolar depression . The genes PLCE1, NOC3L, and CYP2C8 were related to limbic network-FC, all of which have been noted to have a relationship with major depressive disorder (Shi et al., 2012;Q. Wang et al., 2019;Kanders et al., 2020). A previous study investigating the role of PLCE1 in major depressive disorder patients has demonstrated an association with antidepressant remission in female patients, together with other genes within the calcium/calmodulin-dependent protein kinase pathway (Shi et al., 2012). NOC3L eQTLs in the cerebellum and nucleus accumbens have previously been demonstrated to associate with depression severity and antidepressant response (Kanders et al., 2020), and one of the substrates of CYP2C8 is clinically prescribed as treatment for major depressive disorder (phenelzine;Q. Wang et al., 2019). These results seem to suggest that major depressive disorder and antidepressant response involve processes that are impacted by genetic variation in limbic network-FC.
We also find evidence of shared genetic signal in FC across different RSNs. Specifically, we observe a genetically correlated and common genome-wide significant locus for both somatomotor and frontoparietal network-FC near PAX8. PAX8 regulates multiple genes involved in the production of thyroid hormone (Di Magliano et al., 2000), an interesting result considering that both somatomotor and frontoparietal network-FC have been linked to (subclinical) hypothyroidism (Singh et al., 2015;Kumar et al., 2018). Additionally, we detect genetically correlating loci between all RSN-FC and a genome-wide genetic correlation between ventral attention and default mode network-FC. The ventral attention network supports salience processing (Kim, 2010), whereas the default mode network includes areas widespread over the brain and supports emotional processing, self-referential mental activity, and recollection of prior experiences (Raichle, 2015). Increased FC within these two RSNs has been associated with bulimia nervosa (Domakonda et al., 2019) and contributes to episodic memory retrieval (Kim, 2010). Altogether, the shared genetic underpinnings of different RSN-FC that we present here could give a possible explanation how multiple disorders are associated with more than one RSN. We report considerable heritability estimates for RSN-SC (ranging from 10.00% to 15.40%) and identify nine genes that suggest a role for synaptic transmission in the genetics of visual network-SC. For example, STRN encodes for a calmodulin-binding protein that is mostly found in dendritic spines playing a role in Ca 21 signaling (Moqrich et al., 1998), the transcript of FAM175B is a component of a deubiquitylation enzyme complex that has been suggested play a role in synaptic transmission and synaptic plasticity (DiAntonio and Hicke, 2004), and SEMA3A is known as an axonal guidance gene during development (Zhou et al., 2008). The SEMA3A protein has been shown to be upregulated in schizophrenia patients and is suggested to contribute to the developmentally induced impairment of synaptic connectivity in the disorder (Eastwood et al., 2003). Visual network functional hyperconnectivity has been observed in schizophrenia (Damaraju et al., 2014;Ford et al., 2015) and related to visual hallucinations (Ford et al., 2015), but future studies should investigate the equivalent SC component in more detail given our findings.
When investigating the genetic relationship between SC and FC within each RSN, we find no significant genomewide or local genetic correlations. Since the estimation of genetic correlations is dependent on sample size and the heritability estimates of both traits (Bulik-Sullivan et al., 2015), studies with higher power are needed to examine the robustness of these correlation estimates. Future studies could additionally incorporate recent insights that indirect structural connections supporting direct functionally connected regions complicate simple structure to function mapping (Suárez et al., 2020). Our study focused on direct structural connections within RSNs. The possibility that the genetics of RSN-FC overlap with that of indirect pathways that structurally connect brain regions within RSNs via a route beyond the borders of that RSN could therefore be subject to future research. Another possibility is that RSN-FC genetically correlates with RSN-SC if alternatively defined. Here, we have used two metrics to measure properties of connectivity that are most often used in neuroimaging studies and have well established relevance to neuropsychiatry (Buckholtz and Meyer-Lindenberg, 2012;de Lange et al., 2019). These metrics come however with their own limitations (see Kelly et al., 2012;Jones et al., 2013) and alternative metrics that have been used to measure connectivity are the number of streamlines (Jones et al., 2013), index of axonal density (Fieremans et al., 2011), time dependent efficacy of interactions between brain regions (Friston, 2011), or proxy connectivity measures such as morphometric similarity metrics (Seidlitz et al., 2018).
Several limitations must be considered while interpreting our results. First, the definition of RSN used in this univariate GWAS study reduces voxel-level diffusion and functional information to one phenotype by averaging potentially variable connectivity patterns. This could unequally affect more variable higher-order RSNs compared with less variable unimodal RSNs, which would lead to differential statistical power across the RSNs studied here (Helwegen et al., 2023). This could explain why the most significant results are observed in the visual and somatomotor network. Second, it is known that rsfMRI measures are subject to motion distortion, which raises the possibility of differences in measurement error between RSN-FC and RSN-SC. However, given our stringent preprocessing and quality control to enable noise minimization and additional use of rsfMRI-specific covariates in GWAS, we were able to find heritability estimates for RSN-FC that are concordant with previous studies (Roelfs et al., 2022). Third, although UK Biobank provides genetic and uniform MRI data at unprecedented sample sizes, it is evident that even larger sample sizes are needed for discovering the often small genetic effects of polygenic traits (Visscher et al., 2017). The null results observed for some RSN-FC/SC GWAS, partitioned heritability and gene-set analyses might be explained by the multiple comparison correction for the number of phenotypes analyzed, in conjunction with insufficient statistical power. Fourth, some other sample characteristics, such as the European ancestry, age-class and socioeconomic status of subjects, may limit the generalizability of our findings. While we corrected for age and Townsend deprivation index (a proxy of socio-economic status) in our GWAS to reduce this bias, larger and more diverse imaging-genetics datasets are undoubtedly needed.
This study examines the genetic architecture of RSNs, structurally and functionally. We observe several genetic effects for RSNs that highlight relevant biological processes for brain connectivity and related brain disorders. The complexity of structure-function coupling within RSNs is illustrated by the observation that, despite genetic overlap of RSNs within the functional domain, genetic overlap is less apparent within the structural domain and between the functional and structural domains. Altogether, this study advances the understanding of the complex functional organization of the brain and its structural underpinnings from a genetics viewpoint.