Shared genetic influences on resting‐state functional networks of the brain

Abstract The amplitude of activation in brain resting state networks (RSNs), measured with resting‐state functional magnetic resonance imaging, is heritable and genetically correlated across RSNs, indicating pleiotropy. Recent univariate genome‐wide association studies (GWASs) explored the genetic underpinnings of individual variation in RSN activity. Yet univariate genomic analyses do not describe the pleiotropic nature of RSNs. In this study, we used a novel multivariate method called genomic structural equation modeling to model latent factors that capture the shared genomic influence on RSNs and to identify single nucleotide polymorphisms (SNPs) and genes driving this pleiotropy. Using summary statistics from GWAS of 21 RSNs reported in UK Biobank (N = 31,688), the genomic latent factor analysis was first conducted in a discovery sample (N = 21,081), and then tested in an independent sample from the same cohort (N = 10,607). In the discovery sample, we show that the genetic organization of RSNs can be best explained by two distinct but correlated genetic factors that divide multimodal association networks and sensory networks. Eleven of the 17 factor loadings were replicated in the independent sample. With the multivariate GWAS, we found and replicated nine independent SNPs associated with the joint architecture of RSNs. Further, by combining the discovery and replication samples, we discovered additional SNP and gene associations with the two factors of RSN amplitude. We conclude that modeling the genetic effects on brain function in a multivariate way is a powerful approach to learn more about the biological mechanisms involved in brain function.


Abstract
The amplitude of activation in brain resting state networks (RSNs), measured with resting-state functional magnetic resonance imaging, is heritable and genetically correlated across RSNs, indicating pleiotropy. Recent univariate genome-wide association studies (GWASs) explored the genetic underpinnings of individual variation in RSN activity. Yet univariate genomic analyses do not describe the pleiotropic nature of RSNs. In this study, we used a novel multivariate method called genomic structural equation modeling to model latent factors that capture the shared genomic influence on RSNs and to identify single nucleotide polymorphisms (SNPs) and genes driving this pleiotropy. Using summary statistics from GWAS of 21 RSNs reported in UK Biobank (N = 31,688), the genomic latent factor analysis was first conducted in a discovery sample (N = 21,081), and then tested in an independent sample from the same cohort (N = 10,607). In the discovery sample, we show that the genetic organization of RSNs can be best explained by two distinct but correlated genetic factors that divide multimodal association networks and sensory networks. Eleven of the 17 factor loadings were replicated in the independent sample. With the multivariate GWAS, we found and replicated nine independent SNPs associated with the joint architecture of RSNs. Further, by combining the discovery and replication samples, we discovered additional SNP and gene associations with the two factors of RSN amplitude. We conclude that modeling the genetic effects on brain function in a multivariate way is a powerful approach to learn more about the biological mechanisms involved in brain function.

K E Y W O R D S
genetic correlation analysis, genomic SEM, multivariate GWAS, pleiotropy, resting-state networks
However, RSN activity can also be captured by the amplitude of BOLD fluctuations (Bijsterbosch et al., 2017;Zhang et al., 2011), a measure representative of signal changes in the activity of individual RSNs. Previous studies have shown that variations of BOLD amplitude over time in a given brain region are associated with changes over time in its functional connectivity with other regions (Bijsterbosch et al., 2017;Cole, Yang, Murray, Repovš, & Anticevic, 2016). BOLD amplitude-based measures are relevant for explaining human behavior, as demonstrated by their association with cognitive performance (Bijsterbosch et al., 2017;Mennes et al., 2011;Xu et al., 2014;Zou et al., 2013) and life-risk behaviors (Bijsterbosch et al., 2017). More recently, BOLD amplitudes in individual RSNs of the UK Biobank were shown to be (SNP-based) heritable , with estimates on average higher than those scored by (partial) correlation-based measures in the same sample (0:14 < h 2 SNP < 0:36). The genome-wide association study (GWAS) of BOLD amplitude conducted by Elliott et al., 2018 led to the discovery of the first genomic loci associated with individual RSNs: seven RSNs covering prefrontal, parietal and temporal cortices were associated with SNPs in the gene PLCE1; four RSNs mainly covering prefrontal regions were associated with the same three intergenic variants (rs7080018, rs11596664, rs67221163) in chromosome 10; one genome-wide association with a single sensorimotor RSN involved the intronic variant rs60873293.
Next to an overlap in single genetic variants involved in multiple RSNs, significant genetic correlations between different RSNs have been reported using bivariate GWAS analysis  and twin models (Reineberg et al., 2020;Teeuw et al., 2019), with observed correlations between .37 and .79. These results suggest that RSNs are driven by shared genetic variation, indicating the potential for pleiotropy, that is, the same genetic variants being involved in the etiology of different RSNs. One twin study conducted by Reineberg et al., 2020 showed that the heritability of brain connectivity within and across RSNs is represented by three clusters, of which one was defined by low-heritability connections, and two clusters of heritable connections. The latter two can be described as one cluster comprising connections with high heritability in the visual cortex (i.e., "sensory" regions) and a second cluster comprising associations among default mode, frontoparietal, salience, dorsal, and ventral attention regions (i.e., "multimodal association" regions, which integrate inputs from multiple sensory modalities). This broad division of the connectome into sensory networks and multimodal association networks is in line with what was previously found on the basis of clustering of BOLD amplitude across RSNs as well (Bijsterbosch et al., 2017;Zhang et al., 2011). Based on these results, we hypothesize that RSNs genetically diverge according to their "sensory" or "multimodal association" functions.
To identify the SNPs and genes driving this observed pleiotropy multivariate methods can be applied (Grotzinger et al., 2019;Turley et al., 2018;Zhu et al., 2015). Grotzinger et al. (2019) used a novel technique called genomic structural equation modeling (genomic SEM) to model a single genetic factor capturing GWAS associations across multiple psychiatric diagnoses. This multivariate GWAS approach led to the discovery of SNPs that were not observed by any of the separate univariate GWASs of any of the disorders (Grotzinger et al., 2019). In this way, multivariate GWAS provides a new, statistically powerful opportunity to directly characterize the genomic influence on multiple phenotypes simultaneously. Given the observed pleiotropy between brain activation in RSNs, the same approach can be applied to discover the SNPs most associated with shared genetic effects on brain function.
In the current study, we investigated shared genetic etiologies of multiple RSNs within the brain. We used GWAS summary statistics for the amplitude of 21 RSNs throughout the brain made available by the UK Biobank (Bycroft et al., 2018;Elliott et al., 2018;Miller et al., 2016;Sudlow et al., 2015). Our approach was conducted on the GWASs reported for the discovery sample (N = 21,081;Smith et al., 2021), with a replication being carried out on an independent sample from the same cohort (replication sample: N = 10,607 ;Smith et al., 2021). The same approach was then repeated on the GWASs of all the available individuals in both the discovery and replication samples, that is, BIG40 sample (N = 31,688; Smith et al., 2021). First, we estimated the SNP-heritability of the selected RSNs, and for heritable RSNs we modeled their shared genetic structure using genomic SEM (Grotzinger et al., 2019). Next, we performed multivariate GWASs to characterize the SNPs associated with these pleiotropic factors. The multivariate GWAS findings obtained with the BIG40 sample were further interpreted via functional annotation of top GWAS loci and gene-mapping with the functional mapping and annotation (FUMA) tool (Watanabe, Taskesen, van Bochoven, & Posthuma, 2017) and gene-wide and gene-set analysis in MAGMA (de Leeuw, Mooij, Heskes, & Posthuma, 2015). Finally, we also tested whether the newly found genomic factors were genetically correlated with neuropsychiatric and physical traits.

| SNP-based heritability of RSNs
We obtained GWAS summary statistics of BOLD amplitudes of 10 multimodal association and 11 sensory RSNs (Figures 1 and S1, and Tables S1 and S2) measured both in the discovery and the BIG40 samples , with 21,081 and 31,688 adult individuals, respectively (see Section 4.1). With the BIG40 sample, the h 2 SNP of amplitude in all the 21 RSNs was false discovery rate (FDR)-corrected significant (adjusted p [FDR] ≤ .05; Figure S1 and Table S1). Therefore, we kept all 21 RSNs for subsequent analyses on the BIG40 sample.
The heritability estimates ranged between 0.05 and 0.17, with multimodal association networks showing on average higher h 2 SNP than sensory networks (average h 2 SNP = 0.11 and 0.07, respectively). Within the discovery sample, 19 of the 21 RSN amplitudes showed a significant SNP-based heritability, while two sensory networks involved in secondary visual processing (SN2-3) had nonsignificant h 2 SNP estimates (h 2 SNP = 0.036 and 0.038; Figure S1 and Table S2). The two networks were thus excluded from subsequent analyses of the discovery and replication samples.

| Genetic correlations between RSNs
To test the existence of shared genetic etiologies between the heritable RSN amplitudes, we calculated genetic correlations using Linkage Disequilibrium Regression Analysis (Bulik-Sullivan et al., 2015) available within the genomic SEM package (Grotzinger et al., 2019). For more details on the genetic correlation values and respective standard errors and p values, see Table S3. For the genetic correlation results obtained with the discovery sample, consult Table S4.

| Genomic structural equation modeling
To characterize the common underlying genetic etiologies between heritable RSNs, we derived latent genomic factors using genomic SEM (Grotzinger et al., 2019). We chose the most optimal model on the basis of exploratory factor analysis (EFA) in the discovery sample. The EFA results are summarized in Figure 3, showing that the two-factor model explained 55% variance, 16% more variance than the onefactor model, while the addition of a third factor did not explain substantially more variance (Karlsson Linnér et al., 2019;Levey et al., 2020). These results indicate that the most optimal model in representing the pleiotropy among these RSNs consists of two factors. The factor loadings retrieved by this EFA are available in Table S5.
We used confirmatory factor analysis (CFA) to test the model fit of the two-factor model in the discovery sample, and retested the same model in the replication sample (N = 10,607; Smith et al., 2021).
The model fit estimates are reported in Table 1 for the two samples.
F I G U R E 1 SNP-based heritability results obtained for multimodal association and sensory networks. Cortical surface maps displayed show the multimodal association and sensory networks measured in the BIG40 sample. Multimodal association networks are displayed at the top and the sensory networks at the bottom. Both the medial and lateral views of RSNs in the left and right hemispheres are displayed from left to right. RSNs are color-coded according to their SNP-based heritability-proportion of variance in the trait explained by SNP effects-whose scales are displayed on the right. RSN, resting state network; SNP, single nucleotide polymorphism   Table 1). By comparing, in the discovery sample, the chi-square and Akaike information criterion (AIC) statistics between the two sets, we observed that excluding non-significant factor loadings from the model led to lower values retrieved by both statistics, and thus an improved model fit. Further, by testing this model on our replication sample, we observed that the replication sample had an even better model fit compared to the discovery sample. The factor loading results retrieved across the discovery and replication CFA (Tables S6-S8) showed that 11 out of the 17 RSN amplitude associations with the two factors were replicated based on the nominal significance (p < .05) reported with the replication sample (of which five also significant upon Bonferroni correction; p [Bonferroni] ≤ .05/17 = .0029).
The genomic SEM approach conducted on the discovery sample was then repeated on the BIG40 sample. Despite including two additional RSNs compared to the discovery sample (SN2-3), the EFA on the BIG40 sample led to the same optimal number of factors (see where the first factor (F1) comprises all 10 multimodal association networks (MA1-10) and two sensory networks (SN4 and SN10); whereas the second factor (F2) consists of six sensory networks (SN2, SN5-9). These results resemble highly the genomic SEM outcome obtained with the 17 RSN amplitudes included in the discovery sample (Tables S6-S8), which did not include SN2 due to nonsignificant SNP-based heritability. In Tables S9-S11, we include the nominal and Bonferroni-corrected p values of the factor loadings leading to the results in Figure 4. For completeness, we also report results of the one-factor model, which also had a reasonable fit, for each step in the Tables S12-S18.

| Multivariate GWAS results
We estimated the SNP effects driving the pleiotropy of RSNs using multivariate GWASs of the two latent genetic factors. Both F1 and F2 showed significant SNP-based heritability in the discovery sample (F1: showed genome-wide significant associations with F1 (p < 5EÀ8). Table 3 shows the results obtained for the nine independent genomewide significant SNPs of the 142 SNP associations (for GWAS plot, see Figure S2). We found that all nine SNP associations were replicated with their nominal significance p < .05 in the multivariate GWAS on the replication sample (see Table 3). All nine SNP associations would remain replicated if we adopted a more stringent Bonferroni correction accounting for the number of independent genome-wide significant SNPs (p[Bonferroni] = .05/9 = .0056). For F2, no SNP F I G U R E 3 Summary of exploratory factor analysis. Plot displaying the percentage of cumulative explained variance (r 2 ) from up to six-factor models tested using exploratory factor analysis (EFA) on the discovery (left) and BIG40 samples (right) (a); Cumulative explained variance by the one-, two-, and three-factor models tested using EFA on the discovery (left) and BIG40 samples (right) (b); the added explained variance corresponding to an additional factor in the model is shown in parenthesis T A B L E 1 Summary of the two-factor confirmatory factor analysis in the discovery and replication samples Note: For each sample listed in the first column, the second column distinguishes the stages composing the CFA approach. For the discovery sample, CFA consisted of two stages: the first stage, that is, EFA-based model, tested the model design for the two-factor model indicated by our EFA approach; the second stage, that is, modeling for multivariate GWAS, only kept RSNs that showed Bonferroni-corrected factor loadings. The CFA conducted for the replication sample consisted of a single stage, which used the model employed in the multivariate GWAS of the discovery sample. In each stage, a model with a given number of factors was tested (third column), with a given number of RSNs (fourth column). For each model, we display the chi-square statistic, degrees of freedom, AIC, CFI, and SRMR, from the fifth to the ninth columns.
reached genome-wide significance in the discovery sample ( Figure S2), The multivariate GWAS of F1 in the BIG40 sample reported, in addition to the SNPs found with the discovery sample, an additional amount of 357 SNPs, making a total of 498 SNPs circumscribing seven genomic loci (see Table 4 and Figure S3). Of these SNPs, 128 located in three of the seven loci had genome-wide significant Q SNP statistics (p < 5EÀ8), indicating that some SNP effects in these loci are driven by specific RSN, rather than by the multiple RSNs associated with F1. The analysis on the BIG40 sample also revealed that F2 is associated with 21 SNPs with genome-wide significance in a loci with the lead SNP being rs6737318 in chromosome 2 (p = 2.15EÀ9; nearest gene PAX8; Figure S3). However, all these 21 SNPs reported genome-wide significant Q SNP statistics (p < 5EÀ8), and appeared specifically driven by the SNP effects on four sensorimotor networks (SN5-8).
All the genomic loci reported for both factors were also found to be associated (p < 5EÀ8) in at least one of the 18 previous original

| Functional characterization of top GWAS loci
We interpreted our multivariate GWAS results by conducting functional annotation and gene-mapping of genomic loci using FUMA (Watanabe et al., 2017). In addition to the 489 genome-wide significant SNPs reported for F1 with the BIG40 sample, FUMA analysis identified 159 other SNPs in linkage disequilibrium (LD) with these genome-wide significant SNPs, making a total 648 candidate SNPs distributed among seven genomic loci (see Table 4). With the func- SNPs that map to these genomic loci.

| Gene-wide and gene-set results
To investigate whether our multivariate SNP-associations aggregated in a biologically meaningful way, we performed gene-wide and geneset association analyses for F1 and F2 using MAGMA (de Leeuw et al., 2015). Within the BIG40 sample, we found 14 genome-wide significant genes associated with F1 ( Figure (Table S26). For F2, we discovered one gene-wide association for ANO1, but no significant gene-sets (Table S27). Additionally, we investigated via MAGMA tissue expression profile analysis whether the genes associated with F1 and F2 were enriched in 30 general human tissue types and 53 more specific tissue types. No significant enrichment was found for F1 or F2 genes ( Figures S4 and 5).

| Genetic correlations with neuropsychiatric and physical traits
To examine shared genetic effects between the two multivariate RSN factors, estimated using the BIG40 sample, and 10 preselected neuropsychiatric and physical traits, we performed genetic correlation analyses with GWAS summary statistics. The genetic correlation results are reported in Figure 6. Note: For the BIG40 sample, the first column distinguishes the two stages composing the CFA approach, respectively: the first stage, that is, EFA-based model, tested the model design for the two-factor model indicated by our EFA approach; the second stage, that is, modeling for multivariate GWAS, only kept RSNs that showed Bonferroni-corrected factor loadings. In each stage, a model with a given number of factors was tested (second column), with a given number of RSNs (third column). For each model, we display the chi-square statistic, degrees of freedom, AIC, CFI, and SRMR, from the fourth to the eighth columns.
F2 was nominally significantly correlated with genetic factors driving autism spectrum disorder, Alzheimer's disease, body-mass index (BMI), and bone density. For more details on the genetic correlation values and respective standard errors and p values, see Table S28.

| DISCUSSION
We investigated the genomic basis of pleiotropy of brain function in 21 RSNs across the brain. We discovered that two latent genetic factors best captured the genomic influence on the amplitude of RSNs throughout the brain, with 11 RSN amplitude associations replicated in an independent sample. The first factor was associated with multimodal association networks and two sensory networks; the second factor represented only sensory networks. Further, we found that the first factor was associated with SNPs and genes with implications for our understanding of the molecular basis of brain function.
Our genomic factor analyses point to a genetic divergence of multimodal association and sensory functions. This distinction is in line with previous studies using functional connectivity measures (Reineberg et al., 2020) and phenotypic analyses of RSN amplitudes (Bijsterbosch et al., 2017;Zhang et al., 2011). Brain regions involved in sensory and multimodal association functions have also been found to differ in cytoarchitectonic properties (Mesulam, 1998 (Hawrylycz et al., 2012). Together with our findings, the extensive evidence of genetic and brain differences between these two factors may potentially reflect the known differences in the period of maturation between their respective brain regions (Fuhrmann, Knoll, & Blakemore, 2015). In addition, in evolution humans show more pronounced cortical expansion in multimodal association networks than they do in sensory networks compared to other primate species Buckner & Krienen, 2013;Wei et al., 2019).
Thus, differences between sensory and multimodal brain networks have been consistently indicated across biological disciplines from neurodevelopment, to neurophysiology, to evolution. With our findings, we suggest that this divergence between the sensory and multimodal association systems may also be represented by partly distinct effects of common genetic variation in the BOLD amplitude of RSNs.
As more, bigger, and more deeply phenotyped resources looking at genome, brain and their intermediate biology (e.g., epigenome or transcriptome) become available, future studies may test with increased power whether this genomic divergence extends to other levels of biology.
Our first factor is mainly marked by the influence of all the multimodal association RSNs included in the model, covering also the effects coming from two sensory RSNs (SN4 and SN10; Figure 4).
Although SN4 was determined to belong to the sensory system based on the classification carried out by Bijsterbosch et al., 2017, SN4 covers a wide range of brain regions involved in language, that is, language network, and was a priori expected to belong to the multimodal association cluster at the phenotypical level. Therefore, it is not surprising that our first factor is driven by the genetic effects coming from SN4. Less expected was the inclusion of SN10, which is consisted of the supplementary motor area and the striatum. We speculate that the striatum, being closely linked to the frontal cortex and important for executive functions such as working memory (Chudasama & Robbins, 2006;Monchi, Hyun Ko, & Strafella, 2006), explains the high genetic overlap of SN10 with other multimodal association networks, and thus its inclusion in the first factor.
F I G U R E 4 Path diagram of two-factor models reported for the BIG40 sample. Orange circles represent the two latent genetic factors of the two-factor confirmatory factor analysis (CFA F I G U R E 5 Manhattan plot of MAGMA gene analysis findings of latent factors F1 and F2 reported with the BIG40 sample. Gene-wide p values of associations in F1 (top), which comprises genetic effects shared among all 10 multimodal association networks (MA1-10) and two sensory networks (SN4 and SN10); and F2 (bottom), consisted of six sensory networks (SN2 and SN5-9). In each plot, genes located across the 22 autosomes labeled along the x-axis are represented by blue dots, whose position along the y-axis represents the log p value scored by their gene-wide association with each latent factor. The red-dashed horizontal line marks the Bonferroni-corrected significance for the number of genes being tested (p[Bonferroni] ≤ 2.64EÀ6) F I G U R E 6 Genetic correlation matrix comparing the two factors of general brain function with neuropsychiatric and physical health traits. Genetic correlations of two genetic factors (F1 and F2), estimated using the BIG40 sample, with 10 neuropsychiatric and physical traits: attention deficit/hyperactivity disorder (ADHD), autistic spectrum disorder (ASD), bipolar disorder (BIP), major depressive disorder (MDD), schizophrenia, Alzheimer's disease (AD), height, body-mass index (BMI), diastolic blood pressure (DBP), and bone density. Genetic correlations at nominal and false discovey rate (FDR)-corrected significance are respectively labeled with * and ** The first factor of general brain function was associated with a total of 648 candidate SNPs distributed among seven genomic loci.
We also detected 15 gene-wide associations, six of which were not previously detected by the univariate GWAS, and we were able to functionally map 96 additional genes relevant to the study of brain physiology. We thus demonstrate that a multivariate genomic approach has additional value in the search for genetic underpinnings of brain function.
Our gene-wide analysis showed that the first factor was associated with APOE, an important risk-gene for Alzheimer's disease as demonstrated by GWASs Kunkle et al., 2019;Lambert et al., 2013), so as by studies focused on the association between Alzheimer's disease and APOE genotype (Corder et al., 1993;Farrer et al., 1997;Liu, Kanekiyo, Xu, & Bu, 2013). This finding points thus to a possible role of neurodegenerative processes in the first factor. This hypothesis is also supported by other gene-wide associations reported in previous GWASs of Alzheimer's disease, such as MAPK7 (Nazarian, Yashin, & Kulminski, 2019) and APOC1 (Lo et al., 2019;Nazarian et al., 2019), and the functional mapping of LGI1, which was previously reported in relation to beta-amyloid measurement in cerebrospinal fluid (Chung et al., 2018;Li et al., 2015), a biomarker for Alzheimer's disease (Blennow & Hampel, 2003;Frank et al., 2003;Sunderland et al., 2003). An eventual link between this factor and aging-effects potentially reflective of Alzheimer's disease was also suggested by a nominally significant genetic correlation analysis with Alzheimer's disease, opening the possibility that this link is reflected by association patterns at the genome-wide scale.
Interesting significant gene-wide associations also included FHL5, a gene previously associated with migraine (Adewuyi et al., 2020;Gormley et al., 2016), spatial memory (Greenwood et al., 2019), and cerebral blood flow (Ikram et al., 2018), and EPN2, which encodes a protein involved in notch signaling endocytosis pathways, and has previously been associated with educational attainment (Kichaev et al., 2019;Lee et al., 2018) and schizophrenia (Goes et al., 2015); this indicates that notch signaling, known for its role in neurodevelopment and the onset of psychiatric disorders (Hoseth et al., 2018;Lasky & Wu, 2005), may also have an influence on general brain function in adulthood. However, the genes retrieved by our gene-wide analyses were not all related to traits exclusively relevant to the brain, but also to cardiovascular (Ehret et al., 2011;German, Sinsheimer, Klimentidis, Zhou, & Zhou, 2020;Giri et al., 2019;Hoffmann et al., 2017), metabolic (Hübel et al., 2019;Krumsiek et al., 2012;Rask-Andersen, Karlsson, Ek, & Johansson, 2019;Shin et al., 2014), and drug response traits (Cha et al., 2010;Ji et al., 2014;Pardiñas et al., 2019;Takeuchi et al., 2009). BOLD amplitude, being a blood-based measure, may also be susceptible to genetic effects affecting blood-related traits that are not necessarily specific to the brain. The gene-wide result for PLCE1 is an example of such an observation, since it was previously reported for 38 other phenotypes, covering brain (e.g., migraine), cardiovascular (e.g., hypertension, blood pressure), and more general metabolic traits (e.g., BMI). The associa- Despite known associations of the above rfMRI-associated genes with neuropsychiatric and physical traits and previously reported phenotypic associations between these traits and rfMRI-derived imaging phenotypes (Badhwar et al., 2017;Cortese et al., 2021;Lau et al., 2019;Miller et al., 2016;Mulders et al., 2015;Wojtalik et al., 2017), we did not detect genetic correlations between our two genetic factors for brain function and those other neuropsychiatric and physical traits that remained significant after correcting for multiple comparisons. However, the nominal significance reported in six genetic correlations involving neuropsychiatric disorders (major depressive disorder, autistic spectrum disorder, and Alzheimer's disease) and physical traits (BMI and bone density) still suggests that eventual links between general brain function and these phenotypes may be explained by additive effects of common variants across the whole genome.
The discovery of SNPs and genes mainly associated with traits not specific to the brain suggests that other potential sources of genetic signals may drive our multivariate GWAS results of BOLD amplitude. There is a chance that these sources may include typical MRI confounders previously shown to be associated with BOLD amplitude-based measures, such as head motion (Bijsterbosch et al., 2017) and physiological fluctuations associated with respiration or heart functions (Golestani, Wei, & Chen, 2016;Kannurpatti & Biswal, 2008). The association between BOLD amplitude and head motion was particularly noted in sensory RSNs, and it was influenced by the (decreasing) arousal of participants during the MRI scanning (Bijsterbosch et al., 2017), for example, participants tend to become increasingly sleepy as scanning duration increases. The arousal of participants reflects their daily sleep duration and quality, which are also associated with changes in physiological fluctuations measured in cardiac and breathing rate (Snyder, Hobson, Morrison, & Goldfrank, 1964). Taking into account the quality control (QC) implemented in the imaging data leading to the GWAS summary statistics included in our analysis (Alfaro-Almagro et al., 2021), we expect that the confounding effects introduced by both physiological structured noise and head motion are minimized and do not explain our results. Yet, we do not exclude the possibility that residual effects from these variables are still present, which does not necessarily imply the presence of noise in our genomic factors. For example, head motion has been increasingly perceived as a complex measure that also carries behaviorally relevant effects that are heritable (Couvy-Duchesne et al., 2014;Hodgson et al., 2017) and genetically correlated with demographic and behavioral traits (Hodgson et al., 2017). The more recent awareness of arousal as an MRI confounder (Bijsterbosch et al., 2017) leads to a similar scenario, given its relationship with heri-  (Lane et al., 2017;Madrid-Valero, Rubio-Aparicio, Gregory, Sánchez-Meca, & Ordoñana, 2020). The case of arousal applies particularly to our second factor (SN2, SN5-SN9), not only due to the previously reported phenotypic association between sensory RSN amplitudes and arousal (Bijsterbosch et al., 2017), but also because the lead SNP associated with our second factor (rs6737318) is also associated with sleep duration (Dashti et al., 2019; consult Table S25).
This study should be viewed in light of several strengths and limitations. Strengths of our study are the use of GWAS results of large resting-state fMRI samples, which provided the power necessary to run this analysis. Furthermore, we used state of the art novel technologies to find shared genetic etiologies in summary statistics including genomic SEM, which provided statistical power-boosting through the joint analysis of GWASs. Our results provide a new, data-driven basis for studying biological pathways relevant to brain function, by integrating multiple data sources spanning genomics, epigenomics, and transcriptomics. However, this characterization was limited by the data sources that are currently available. As more resources become publicly available and integrated in FUMA and equivalent platforms, in the future an even broader genetic mapping of traits will be possible.
Another limitation of our study is the fact that our approach focused exclusively on the effects of common SNPs, without including the effects of rare genetic variants or gene-environment interactions and correlations. Including rare variation in follow-up studies and more extended explicit modeling of gene-environment interplay may provide even more insight into the biological pathways underlying brain function.
In conclusion, we show that pleiotropy in heritable RSNs is best represented by a two-factor model mainly distinguishing the genetic influences on multimodal association from those on sensory networks.
GWAS-based analysis of these genetic factors led to the discovery of relevant SNPs and genes. With our approach, we demonstrate that taking advantage of the pleiotropy of RSNs using multivariate genome-wide approach provides new insights in the genetic and molecular roots of brain function.

| GWAS sample
We used GWAS summary statistics from the UK Biobank initiative, publicly available in a second release via Oxford Brain Imaging Genet-

| Description of RSNs
The 21 RSNs covering spontaneous BOLD fluctuations in the brain were labeled based on the clustering analyses conducted in Bijsterbosch et al., 2017, which appointed RSNs to one of two distinct system categories: multimodal association and sensory systems. In Table S29, we show the system category given to each RSN, the respective UK Biobank label, and its respective two-dimensional anatomical visualization. Visualization of RSNs is also provided by the UK Biobank online resources (fmrib.ox.ac.uk/ukbiobank/group_means/ rfMRI_ICA_d25_good_nodes.html). In the following step, the covariance matrices estimating the pleiotropy among heritable RSN amplitudes were retrieved using the multivariate extension of LDSC distributed by the genomic SEM package.

| Genomic structural equation modeling
We obtained (a) a genetic covariance matrix quantifying the genetic overlap among the RSNs; (b) the respective matrix containing the standardized genetic covariance values (i.e., genetic correlations); (c) a sampling covariance matrix informative of the standard errors associated with the genetic covariance measures.
To determine the number of factors in the model, and which imaging phenotype loaded on which factor, we conducted an EFA with maximum likelihood estimation. Before running EFA, the LDSCderived covariance matrix was smoothed to the nearest positive, as part of the default genomic SEM pipeline. We tested EFA with one factor and repeated the same step for an increasing number of factors up to six. We selected the highest number of factors leading to an explained variance increase (r 2 ) of equal or more than 10% (Levey et al., 2020). For all the modeling results, positive or negative factor loadings with magnitudes equal or higher than 0.35 were assigned to a given factor, identical to Grotzinger et al. (2019).
For the most optimal model, we ran CFA using the genomic SEM package, in order to estimate the factor loadings of the variables included in the model and evaluate the respective model fit. Both the genetic and sampling covariance matrices were analyzed using weighted least squares estimation, providing fit statistics and inferred factor loadings. We retained factor loadings at a Bonferroni significance level across the factor loadings within the model (p [Bonferroni] ≤ .05/number of factor loadings). Further, with the model retaining Bonferroni-significant factor loadings in the discovery sample, we conducted a CFA on the independent replication sample, in which the replication of the factor loadings was determined by "nominal" significance (p ≤ .05). reporting genome-wide associations with SNPs located in that same locus. Gene-mapping was performed by (a) selecting genes located within 10 kb of each SNP, (b) annotating SNPs based on their eQTL enrichment in the data resources listed in Table S30, and (c) the chromatin interactions depicted in the HI-C data resources reported in

| Gene-wide and gene-set analyses
To test for aggregated association of multiple SNPs within genes, we performed gene-wide analyses on the multivariate GWAS results of the BIG40 sample. We then performed gene-set analysis for curated genesets and GO terms from MsigDB c2 and c5, respectively, testing for the presence of pathways associated with these factors. Furthermore, we performed tissue gene expression analysis in the genomic factors. These analyses were all performed using the MAGMA v1.08 software (de Leeuw et al., 2015) as embedded within the FUMA platform (Watanabe et al., 2017), for details see the Supplementary Methods.

| Genetic correlations with other traits
To examine shared genetic effects between the RSN genomic factors estimated with the BIG40 sample and other traits, we performed genetic correlation analyses with GWAS summary statistics from 10 selected traits. We followed the same QC and bivariate genetic analysis procedures used for RSN amplitudes (see Section 4.3). Among the selected GWAS summary statistics, we included six neuropsychiatric disorders with high prevalence in the population that are widely associated to RSN function in literature (Badhwar et al., 2017;Cortese et al., 2021;Lau et al., 2019;Mulders et al., 2015;Wojtalik et al., 2017).

ACKNOWLEDGMENTS
The research leading to these results received funding from the Bralten gratefully acknowledges funding from the NWO Innovation program (Veni 09150161910091). We would like to thank Andrew D.
Grotzinger (University of Texas at Austin) for the support provided in conducting the genomic SEM approach.

CONFLICT OF INTEREST
B. Franke has received educational speaking fees from Medice. C. F.
Beckmann is director and shareholder of SBGneuro Ltd.

AUTHOR CONTRIBUTIONS
The conception of the idea motivating this study was elaborated by João Guimaraes, E. Sprooten, and J. Bralten. João Guimaraes performed the analyses, which were supervised by J. Bralten and E. Sprooten. João Guimaraes, E. Sprooten, and J. Bralten wrote the paper with contributions from the remaining coauthors.

DATA AVAILABILITY STATEMENT
The GWAS summary statistics on the amplitude of 21 RSNs are publicly available in a second release from the UK Biobank initiative via Oxford