Genomic Structural Equation Modeling Reveals Latent Phenotypes in the Human Cortex with Distinct Genetic Architecture

Abstract Genetic contributions to human cortical structure manifest pervasive pleiotropy. This pleiotropy may be harnessed to identify unique genetically-informed parcellations of the cortex that are neurobiologically distinct from functional, cytoarchitectural, or other cortical parcellation schemes. We investigated genetic pleiotropy by applying genomic structural equation modeling (SEM) to map the genetic architecture of cortical surface area (SA) and cortical thickness (CT) for the 34 brain regions recently reported in the ENIGMA cortical GWAS. Genomic SEM uses the empirical genetic covariance estimated from GWAS summary statistics with LD score regression (LDSC) to discover factors underlying genetic covariance, which we are denoting genetically informed brain networks (GIBNs). Genomic SEM can fit a multivariate GWAS from summary statistics for each of the GIBNs, which can subsequently be used for LD score regression (LDSC). We found the best-fitting model of cortical SA identified 6 GIBNs and CT identified 4 GIBNs. The multivariate GWASs of these GIBNs identified 74 genome-wide significant (GWS) loci (p<5×10 -8 ), including many previously implicated in neuroimaging phenotypes, behavioral traits, and psychiatric conditions. LDSC of GIBN GWASs found that SA-derived GIBNs had a positive genetic correlation with bipolar disorder (BPD), and cannabis use disorder, indicating genetic predisposition to a larger SA in the specific GIBN is associated with greater genetic risk of these disorders. A negative genetic correlation was observed with attention deficit hyperactivity disorder (ADHD), major depressive disorder (MDD), and insomnia, indicating genetic predisposition to a larger SA in the specific GIBN is associated with lower genetic risk of these disorders. CT GIBNs displayed a negative genetic correlation with alcohol dependence. Jointly modeling the genetic architecture of complex traits and investigating multivariate genetic links across phenotypes offers a new vantage point for mapping the cortex into genetically informed networks.


INTRODUCTION
A number of different neurobiological markers have been employed in conjunction with various organizational schemes to map the human cortex.It is possible that individual differences in regional cortical surface area (SA) and cortical thickness (CT) may drive factors that affect each person and each region independently.However, the covariance structure of regional SA and CT reveals that individual differences are systematically coordinated within communities of brain regions, uctuate in magnitude together within a population, may be instantiated as structural covariance networks (SCN) 1 , and partially recapitulate established organizational schemes [2][3][4][5] .For instance, SCN organization is consistent with topological patterns of cortical maturation observed throughout developmental stages from childhood and adolescence into early adulthood 6 , and the same patterns are then targeted by neurodegenerative diseases in late life 7,8 .Second, brain regions with highly correlated CT or SA often represent networks that perform dedicated cognitive processes 1,9,10 .Third, regions within SCNs tend to be directly connected by white matter tracts.Indeed, about 40% of SCN connections show convergent white matter ber connections, although other relationships captured by SCNs are distinct from ber connectivity 5 .
The correlation structure between regions represented by a SCN is in uenced by both the environment and genetics.The genetic factors underlying structural correlations closely resemble functional and developmental patterns 4,5,11 .We will refer to these patterns of genetic correlations between brain regions as genetically informed brain networks (GIBNs).Genetic correlations of CT and SA regions have been examined in twin studies 12,13 .These genetic in uences were recapitulated in over 400 twin pairs, to show that the cortex is organized genetically into communities of structural and functional regions, is hierarchical, modular, and bilaterally symmetric 11 .Their genetically informed parcellation identi ed 12 spatially contiguous regions that qualify as GIBNs.
While twin studies have laid important groundwork regarding genetic correlations of the brain, they have several limitations.First, twin studies do not provide speci c genetic variants associated with genetically correlated regions 11 and therefore offer an incomplete characterization of cortical pleiotropy.Second, twin studies rely on the equal environment assumption, which may be invalid for some traits.Third, quantifying the genetic correlation between CT/SA and assembling a well-powered cohort of low prevalence traits such as schizophrenia (0.5% prevalence) 14 or bipolar disorder (1% prevalence) 15 is extremely di cult due to the rarity of pairs where twins affected by one or both traits.Recently, genetic correlations between brain regions derived from GWAS results have been applied to estimate the contribution of common genetic variation to CT/SA heritability 16 .This method confers several advantages over twin studies as they do not have the same assumptions, allow effect-size estimation for individual variants, and have the ability to test genetic correlations with other traits in different populations.These SA and CT GWAS results reveal pleiotropy and genetic correlation across many neuroimaging phenotypes 17,18 .
Genomic structural equation modeling (gSEM) is a multivariate statistical method that can leverage the genetic architecture of multiple genetically correlated phenotypes to derive relatively few latent phenotypes, which describe the observed genetic correlation and elucidate loadings of multiple phenotypes onto the latent phenotype 19 .Therefore, gSEM applied to GWAS offers a genetically informed clustering of the cortex that may be neurobiologically distinct from functional and cytoarchitectural parcellations 6,20 .Multiple regions that have signi cant loadings on a particular factor de ne the brain regions that can be described as a GIBN.Importantly, gSEM can be used to estimate the strength of association between genetic variants and each latent factor in a multivariate GWAS of each GIBN using GWAS summary statistics for the individual traits.Thus, gSEM provides a description of the underlying genetic architecture of the traits being examined and effect size estimates for speci c SNPs and their association with the latent factors.
In the present study, we sought to elucidate the genetic architecture of 34 regional SA and CT phenotypes reported in the ENIGMA-3 GWAS of over 50,000 primarily healthy individuals.We hypothesized that gSEM might identify cortical SA networks consistent with the 12 clusters described by Chen et al. 11 , along with other viable solutions.The genetic correlations reported in Grasby et al. 18 , were stronger within major anatomical lobes than across lobes.Thus, while we predicted gross lobar structure may be re ected by GIBNs, we further predicted that GIBNs would re ect the complex relationships captured by functional networks, canonical resting-state networks (RSN), ber tract networks, and other neurobiological systems 6,11 .We hypothesized that most genetic variants discovered by the ENIGMA-3 cortical GWAS would in uence GIBNs.We also sought to discover novel genetic markers and links between known genetic variants and GIBNs.
Our motivation for the present analysis was that there is robust evidence of disrupted cortical structure and function for most psychiatric disorders 21,22 .We also know psychiatric disorders are polygenic and there is signi cant genetic correlation between disorders [23][24][25][26][27] .We hypothesize that genetics may in uence psychiatric illness by initially acting to disrupt brain structures, or impart vulnerabilities that lead to disrupted brain structures in the presence (or absence) of speci c environmental exposures [28][29][30] .We further hypothesized genetic correlations between GIBNs and major neuropsychiatric disorders that are stronger than correlations between global measures of SA or CT.

Data
We used the results of the ENIGMA-3 cortical GWAS as reported in Grasby et al. 18 that identi ed genetic loci associated with variation in cortical SA and CT measures in 51,665 individuals primarily (~ 94%) of European descent, from 60 international cohorts.In that study, phenotype measures were extracted from structural MRI scans for 34 regions de ned by the Desikan-Killiany atlas using gyral anatomy, which establishes coarse partitions of the cortex 31 .This study analyzed global measures of total cortical SA and average CT, as well as 34 regional measures of SA and CT averaged across left and right hemisphere structures to yield 70 distinct phenotypes.Multiple testing correction in the ENIGMA-3 GWAS was based on 70 independent phenotypes with a GWS threshold of P ≤ 8.3×10 − 10 .We accessed the GWAS summary results for the 34-regional bilateral analyses.The primary GWAS analyses presented in Grasby et al. adjusted for global SA and mean CT.However, we utilized alternate results without global adjustments, as the global-adjusted results produce multiple negative genetic correlations between regions, which might be artifactual and lead to uninterpretable factor loadings.Regional SA and CT metrics were analyzed separately due to computational limitations and because negative genetic correlations between SA and CT could complicate model interpretation 16,18,32 .

Analysis
Our analyses were performed using the genomic-SEM R package 19 .The gSEM was performed twice, once for 34 SA regions and once for 34 CT regions.The gSEM tting process includes an exploratory factor analysis (EFA) stage and a con rmatory factor analysis (CFA) stage.To avoid over tting, we analyzed odd chromosomes in the EFA and even chromosomes in the CFA.Whereas SEM often ts multiple models corresponding to a priori hypotheses built on theoretical models, we took a hypothesis free (data driven) approach.In the EFA, we t models allowing for 1 to 10 factors, for each of SA and CT.Scree plots were examined to ensure that 10 factors would be su cient (see Figures S1 and S2).In the EFA step, positive factor loading estimates greater than a pre-speci ed threshold from the EFA were carried forward to the CFA to be re-estimated, and the remaining loading parameters were set to zero 33 .As there was no consensus on factor loading cutoffs 19,34 , we tested two thresholds: 0.3 and 0.5.Cross loadings were allowed if they exceeded the threshold.Factors that loaded on only a single region were removed as single regions do not constitute factors.Therefore, some models with a large number of factors ended up as redundant and were not carried forward to CFA as the investigation of single regions was already carried out by Grasby et al 18 The Akaike Information Criteria (AIC) was used as our primary measure of model t.For our purposes, a model which minimized the AIC was deemed optimal.Standardized root-mean square residual (SRMR), model 2 , and Comparative Fit Index (CFI) were also calculated.As opposed to regression modeling, where signi cant statistics represent the strength of association between the predictors in the response, in SEM modeling, a signi cant 2 statistic represents a lack of t.However, with large sample sizes, the 2 statistics can be signi cant regardless of the model, which is not informative.We found all 2 statistics were highly signi cant (p ~ 0), and therefore not reported.
The top-performing factor models in the CFA were further optimized by successive removal of nonsigni cant factor loadings, which is considered standard practice 35 .We additionally t a bifactor model as part of the CFA step to account for the observed correlation between the factors.Speci cally, we t a bifactor model where a "total" CT or SA factor was added, which loaded on all regions and a multi-level model where all EFA factors loaded onto a 2nd order factor.The bifactor models failed to converge and the multilevel models failed to improve model t in all cases; hence these results are not reported.

GIBN Overlap with Alternate Networks and parcellations
To explore the possible relevance of GIBNs to other parcellations of the cortex, we used Dice's Coe cient to measure percent volume overlap.We used permutation testing to determine the signi cance for each Dice's coe cient by estimating the probability that the magnitude of overlap occurred by chance.We used 1,000 iterations of populating a given network with randomly selected brain regions, calculating its Dice's coe cient relative to the parcellation of interest, and then comparing the GIBNs true Dice's coe cient to the null distribution of 1,000 Dice's coe cients.The relative position of the Dice's coe cient for a particular GIBN-to-parcellation comparison within the probability distribution provided the signi cance.False Discovery Rate (FDR) was used to correct for multiple testing with GIBNs, receptors, networks, and clusters.

Multivariate GWAS Analysis
Using the GIBNs from our best tting model, we used gSEM to generate a multivariate GWAS of each GIBN.The GWS associations (p < 5x10 − 8 ) for each GIBN were compared to the signi cant SNPs reported by Grasby et al. with and without the global correction.The FUnctional Mapping and Annotations (FUMA) package 38 was used to annotate results from each GIBN GWAS, including annotating SNPs to speci c genes, identifying independent loci, and identifying potential functional variants.FUMA was run using LD in the 1000G Phase3 EUR reference panel and the default FUMA parameters.
While CT and SA were examined separately, both for computational limitations and conceptual reasons, we used the multivariate GWAS results to estimate the genetic correlation between CT GIBNs and SA GIBNs, hypothesizing that they would be consistent with the negative genetic correlation between average CT and total SA 16 .Additionally, to examine the degree to which the CT and SA GIBNs genetically resembled the overall CT and SA measures, we estimated the genetic correlation between each GIBN and the average CT and total SA (uncorrected for ICV) as reported in the Grasby et al.

Polygenicity Analysis
We examined the signi cant SNPs from the GIBN GWAS, as well as SNPs in LD using FUMA to test for functional associations with established behavioral traits and major neuropsychiatric disorders.First, we examined whether observed variants from the GWAS recapitulated GWS SNPs from previous GWASs of neuroimaging traits including cortical GWASs and other structural neuroimaging parameters 17,39−44 .We also looked for SNPs that were signi cant in GWASs of 12 neuropsychiatric disorders from the Psychiatric Genomics Consortium (PGC): ADHD 45 , alcohol dependence 46 , anorexia nervosa 47 , autism spectrum disorder 48 , bipolar 49 , cannabis use 50 , MDD 51 , obsessive compulsive disorder (OCD) 52 posttraumatic stress disorder (PTSD) 53 , schizophrenia 54 , Tourette's syndrome 55 , and anxiety 56 .Finally, FUMA was used to functionally annotate loci that overlapped with previously published GWAS results.

Genetic Correlation with Psychopathology
We used cross-trait LDSC to identify links between psychiatric disorders and CT-derived GIBNs as well as psychiatric disorders and SA-derived GIBNs 57 .We estimated the genetic correlation between CT-and SAderived GIBNs and neuropsychiatric disorders using their GWAS summary statistics 57 .To limit our need for a multiple testing correction, we limited our analyses to the 12 neuropsychiatric disorders noted above.A false discovery rate corrected p-value (P FDR ) was used to correct for the number of GIBNs (10)   and disorders (12).

Data and Code Availability
GWAS summary statistics used in this paper are available on the ENIGMA consortium website (http://enigma.ini.usc.edu/research/download-enigma-gwas-results).The Genomic SEM package used to analyze the data is publicly available at https://github.com/GenomicSEM/GenomicSEM.The ldsc package is publicly available at https://github.com/bulik/ldsc.The results of the multivariate GWASs of the CT-and SA-derived GIBNs are available at https://pgc-ptsd.com/ about/workgroups/imagingworkgroup/.

Model Fit
The SA-derived 6-GIBN solution resulted in the best overall model t to the genetic covariances generated from the GWAS summary statistics (AIC = 22,712,604, CFI = 0.920, SRMR = 0.062).See Supplementary Table S1 for t statistics for each evaluated model.The 6 SA-derived GIBNs (SA1-SA6) loaded on 24 of the 34 brain regions 18 .The standardized estimates for the 6 SA-derived GIBN models (standardized factor loadings) are presented in Supplementary Table S2 and Fig. 1a.The GIBNs generally encompass contiguous brain regions and many correspond to known neuroanatomical features.SA1 contains loadings for inferior temporal, isthmus cingulate, postcentral, precuneus, superior parietal, supramarginal, and temporal pole.SA2 contains loadings for caudal anterior cingulate, caudal middle frontal, medial orbitofrontal, paracentral, and rostral anterior cingulate.SA3 contains loadings for banks superior temporal sulcus (STS), inferior parietal, and middle temporal.SA4 contains loadings for pars opercularis, pars orbitalis, and pars triangularis, SA5 contains loadings for cuneus, lateral occipital, lingual, and pericalcarine, and SA6 corresponds to the auditory cortex.The 6-factor model indicated substantial correlation between GIBNs (r g =0.61 to 0.91) as reported in Supplementary Table S3.
Factor diagrams for SA-and CT-derived GIBNs are presented in Fig. 2. Consistent with prior work, the SAderived GIBNs were largely distinct from CT-derived GIBNs, although some regional overlap exists.For example, SA5 and CT4 are both 4-region GIBNs, with 3 overlapping regions.

GIBN Overlap with High/Low Neuroreceptor Density Regions
We examined the overlap between CT and SA GIBNs and regions of highest (top 20%; Fig. 5a) and lowest (bottom 20%, Fig. 5b) neuroreceptor densities.We found that CT1 overlapped a region of high neuroreceptor densities for many types of neuroreceptors and a region of low fDOPA receptor density.CT2 and SA2 overlapped regions of high 5-HT1a, 5-HT4, and 5HTT receptor density.SA5 overlapped the high 5HTT receptor-density region.(Fig. 5a).

GWAS of GIBNs
To identify speci c genetic variants that may be in uencing the GIBNs, we performed a multivariate GWAS on each SA-and CT-derived GIBN.Manhattan plots for SA-and CT-derived GIBN GWASs, their associated quantile-quantile (QQ) plots, and genomic in ation factors (λ) are provided in Figures S3 to S12.We observed moderate p-value in ation (λ values between 1.06 and 1.16).However, the single-trait LD Score regression intercepts for SA-and CT-derived GIBNs were all less than 1.02, indicating that the apparent in ation was likely due to polygenicity.A total of 5,843 GWS (p < 5x10 − 8 ) variants were associated with the GIBNs.FUMA 38 mapped these variants to 74 independent regions, including 64 loci associated with the 6 SA-derived GIBNs and 10 loci associated with the 4 CT-derived GIBNs.A phenogram 58 of the genetic associations is presented in Fig. 6.A list of all GWS loci is provided in Table S7.Except for two novel SNPs, all others were previously identi ed in Grasby et al. 18 in either the analyses adjusted for global SA/CT or the unadjusted analyses.The rst novel SNP, rs3006933, near the genes SDCCAG8 and AKT3 on chromosome 1, was associated with SA1 (p = 4.08x10 − 9 ).The other novel SNP, rs1004763, on chromosome 22 in the vicinity of the gene SLC16A8, was associated with CT2 (p = 3.41x10 − 08 ).Notably, many of the 75 GWS loci associated with GIBNs were not associated with global measures, but only with individual CT/SA regional measures, and half (37 out of 75) were more signi cantly associated with the GIBNs than the corresponding global measures.Using FUMA we found no signi cant enrichment of a particular tissue type in either CT-or SA-derived GIBNs and no enriched expression of developmental genes or regulators.

Genetic Correlation between CT and SA
Although CT and SA regions were analyzed separately, we examined the genetic correlation between CT and SA using LDSC analysis of the GIBN GWAS results.The mean genetic correlation between SA GIBNs and CT GIBNs is -0.22 (-0.43 to -0.08; Table S10), whereas the mean genetic correlation between the 6 SA GIBNs is 0.77 (0.61 to 0.91; Table S3) and the mean genetic correlation between the 4 CT GIBNs is 0.76 (0.71 to 0.87; Table S6).The dramatically lower correlation between CT and SA compared to within SA and compared to within CT GIBNs supports separate gSEM analyses of CT and SA phenotypes.

LDSC analysis of Genetic Correlation
We examined the genetic correlation between CT and SA GIBNs and psychiatric disorders.The LDSC analysis of SA GIBNs is reported in Table S8.ADHD exhibited a signi cant negative genetic correlation with all SA-derived GIBNs except SA4 (r g =-0.13 to -0.20,p = 3.29x10 − 6 to 0.0038,p FDR =0.00040 to 0.039).
The global measures for thickness and SA are genetically correlated to many of the same psychopathology traits as the GIBNs, but the genetic correlations with global measures are almost always less signi cant than the genetic correlations with the most strongly associated GIBNs.The details of these results are provided in Tables S8 and S9.

DISCUSSION
The goal of the present study was to leverage the pleiotropic architecture of the human cortex to investigate genetic factors underlying CT and SA, and to identify further links between the genetics of CT, SA, and psychopathology.We applying gSEM to jointly model the genetic architecture of 34 brain regions using results from the ENIGMA-3 GWAS 18 .The process was undertaken with gSEM to generate several possible solutions, from which the best-model t was selected.This solution organized brain regions to optimally assign genetic pleiotropy to 6 SA-and 4 CT-derived latent factors, which we have termed genetically informed brain networks (GIBNs).
The GIBNs we generated may be compared to similar structures generated from twin studies.Using 400 twin pairs, Chen et al. generated twelve genetically-informed clusters from vertex-based SA measures 11 .
Chen et al. 11 reported heritability estimates and genetic correlations between genetically informed parcels that are more consistent with classical anatomically-de ned sulcal and gyral boundaries, Brodmann de nitions, and cytoarchitectural patterns than our GIBNs.However, in one respect, the present results are more informative than Chen et al. 11 , as we provide SNPs associated with each of the GIBNs.
The overlap of GIBN GWS loci with prior GWAS of neuroimaging phenotypes or psychiatric disorders rmly points to the relevance of GIBN-related variants to brain structure and cognition.First, we note that novel variant rs3006933 has been previously associated with subcortical volumes 59 .Novel variants rs3006933 and rs1004763 17,42 have been associated with neuroimaging phenotypes of corpus callosum white matter microstructure 60 .A comparison of our GIBN GWAS with published psychiatric disorder GWAS results found that multiple SNPs linked to SA-derived GIBNs were also implicated in a GWAS of schizophrenia 54 .Speci cally, we identi ed a cluster of 4 loci in the CRHR1 gene strongly associated with SA-derived GIBNs (rs62057153 associated with SA1) in our GWAS (p = 5.22x10 − 17 to 8.45x10 − 21 ).We also observed an association between CT1 and rs11692435 (p = 1.17x10 − 12 ), a schizophrenia-related locus, within the ACTR1B gene.Finally, CT-and SA-derived GIBNs were associated with schizophrenia risk variants in the SLC39A8 gene; namely rs13107325 was associated with CT5 and rs13135092 was associated with SA5.No other traits had GWS variants associated with any of the GIBNs.
Many GIBN-associated SNPs have been associated with other cognitive, behavioral, neuroanatomical, neurofunctional, and neuropsychiatric phenotypes.In addition to rs3006933 noted above 16 , SA6-linked locus rs9909861 61 and SA5-linked SNP rs7570830 16 have been associated with subcortical volumes.Multiple loci associated with SA-derived GIBNs that encompass temporal, parietal, and temporo-parietal association cortices, including the SA1-linked locus rs10109434 62 , the SA3-linked SNP rs2299148 63 , and the SA6-linked locus rs9909861 63-67 have been implicated in academic attainment and cognitive ability.
Resting-state functional connectivity evinces robust patterns of synchronous activity that intrinsically organize into canonical networks 84 .We found that SA-derived GIBNs overlap with several canonical RSNs, such as visual network and SA5 (Dice's coe cient = 0.424), which is composed of cuneus, lateral occipital, lingual, and pericalcarine cortices 36 .Twin-based non-linear multidimensional heritability estimates are among the highest for the visual network (left h 2 m = 0.53; right h 2 m = 0.45) and auditory network (left h 2 m = 0.44; right h 2 m = 0.60) 85 .SA6, which includes superior and transverse temporal cortices, overlaps the auditory cortex from twin-derived genetic parcels (Dice's coe cient = 0.211; Fig. 5).
The functional specializations of the human auditory cortex 36 , which include parts of the lateral prefrontal cortex, Broca's area, and subcentral regions, are needed for human vocalization and language 86,87 .The dorsal attention network (DAN), which directs voluntary allocation of attention, has substantial overlap with SA1 (Dice's coe cient = 0.232, Fig. 5) that is comprised of superior parietal, supramarginal, postcentral, precuneus, isthmus cingulate, and inferior temporal regions.A noteworthy omission from SA1, an important feature of the DAN, are the frontal eye elds (FEF) 88 .Since FEF is not a FreeSurfer parcellation output, it may be poorly represented in ENIGMA cortical GWAS.The DAN has relatively high twin heritability estimates (left h 2 = 0.45; right h 2 = 0.40) 85 .SA4 partially overlaps (Dice's coe cient = 0.259) the frontoparietal network (FPN), which includes pars opercularis, pars orbitalis, and pars triangularis, but lacks the critical temporoparietal structures 89 (Fig. 6).
Behavioral traits and neuropsychiatric disorders showed distinct genetic correlations with SA-derived GIBNs that differ markedly from correlations with CT-derived GIBNs.CT3, located in the middle and superior temporal cortices, and CT4, located in the visual perceptual cortex, were strongly negatively correlated with alcohol use disorder.This divergent relationship between CT-derived and SA-derived networks is consistent with the ENIGMA-3 cortical GWAS where a similar pattern of positive and negative correlations between total brain SA and behavioral traits/disorders was found, but average CT correlations with behavioral traits/disorders were non-signi cant 18 .Speci cally, the ENIGMA-3 GWAS found that total SA was signi cantly positively correlated with cognitive function, educational attainment, Parkinson's disease, and anorexia nervosa, but signi cantly negatively correlated with MDD, ADHD, depressive symptoms, neuroticism, and insomnia.In addition, the SA-derived GIBNs showed distinct genetic relationships to several psychiatric disorders.Several SA-derived GIBNs (SA1, SA2, SA4, SA5) were positively correlated with bipolar disorder, whereas SA-derived GIBNs (SA1, SA2, SA3, SA5, SA6) were negatively correlated with MDD, buttressing prior evidence that MDD and Bipolar are distinct conditions with diverging genetics 27 .While the relationship between these SA-derived GIBNs and MDD converge with the ndings from the ENIGMA total SA results, the relationship with bipolar disorder was novel.Thus, GIBNs may provide additional power to detect genetic relationships when their strength across cortical regions is heterogenous.
Interestingly, although several GIBN-associated SNPs were associated with schizophrenia, no GIBNs were signi cantly genetically correlated with schizophrenia (r g =0.029 to 0.034; p-values > 0.30).While this may be counterintuitive, genetic correlation between phenotypes predicts an overlap in SNPs, but the reverse may not be true.A genetic correlation could be zero when many variants affect both traits, but the direction of effects are uncorrelated across variants 90 .
There is ample evidence that genetic variants that in uence SA are distinct from genetic variants that in uence CT 18 .The results of Panizzon et al 32  By contrast, Panizzon et al 32 focused on genetic correlation from twin data, which means genetic marker associations were not available.Grasby et al 18 was focused primarily on GWAS results and secondarily on reporting genetic correlations between CT and SA.Indeed, the results of van der Meer 16 are completely consistent with the results of Grasby et al 18 with the former reporting a genetic correlation between SA and CT of -0.26 and the latter of -0.32.However, van der Meer et al 16 reports that the 4016 out of 7941 causal variants are shared between CT and SA.Therefore, the relatively high genetic overlap but low genetic correlation is due to a mixture of opposing and agreeing effects from variant.Genetic variation affecting gene regulation in progenitor cell types, present in fetal development, affects adult cortical SA 91 .An increase in proliferative divisions of neural progenitor cells leads to an expanded pool of progenitors, resulting in increased neuronal production and larger cortical SA, which is more prevalent in gyrencephalic species (e.g.humans, primates) 92 .By contrast, loci near genes implicated in cell differentiation, migration, adhesion, and myelination are associated with CT.Our ndings suggest this distinction holds for SA-derived compared to CT-derived GIBNs.We hypothesize that the unique genetic correlations of SA-derived GIBNs and CT-derived GIBNs with behavioral traits/disorders may be explained by the distinct developmental functions of their associated genes 93 .

Limitations
A number of limitations deserve consideration in interpreting the present ndings.The present gSEM was based on the GWAS results of Grasby et al. 17,18 , which averaged left and right hemisphere phenotypic measures.Additionally, Grasby et al 17,18 examined the 34 cortical regions as de ned by the Desikan-Killiany atlas.A high-resolution GWAS of the cortex would allow more exibility and rede ning parcellation boundaries informed by genetic pleiotropy.However, simultaneously analyzing GWAS results for each of 100s or 1000s of vertices would present a computational challenge.Therefore, the GIBNs are a proof-of-concept that genetic correlation can be used to enhance the interpretation of high-dimensional GWAS results and provide novel insights into the relationship between neuroimaging phenotypes and psychiatric disorders.

Conclusion
We harnessed the pervasive pleiotropy of the human cortex to realize a unique genetically-informed parcellation that is neurobiologically distinct from functional, cytoarchitectural, and other established cortical parcellations, yet harbors meaningful topographic similarities to other network schemas.Strong genetic correlation between GIBNs and several major neuropsychiatric conditions, coupled with clear con rmation that nearly all GIBN-associated SNPs play a role in cognitive, behavioral, neuroanatomical, and neurofunctional phenotypes, begins to expose the deeply interconnected architecture of the human cortex.Applying gSEM to model the joint genetic architecture of complex traits and investigate We thank all members of the respective site laboratories within the ENIGMA who contributed to general study organization, recruitment, data collection, and management, as well as subsequent analyses.Most importantly, we thank all of our study participants for their efforts to take part in this study.The funding agencies had no part in the analysis of data or approval of the nal publication.The views expressed in this article are those of the authors and do not necessarily re ect the position or policy of the Department of Veterans Affairs or the US government.
and Grasby et al 18 compared to van der Meer et al 16 are focused on different, albeit related measures.van der Meer et al 16 primarily focused on the overlap of individual genetic variants associated with CT and SA and only secondarily on their genetic correlation.
multivariate genetic links across phenotypes offers a new vantage point for mapping genetically informed cortical networks.Declarations SUPPORT National Institute for Mental Health Grant No. R01-MH111671, R01-MH129832, and VISN6 MIRECC (to RAM); VA Merit Grant Nos.1I01RX000389-01 (to RAM) and 1I01CX000748-01A1 (to RAM); National Institute of Neurological Disorders and Stroke Grant Nos.R01-NS086885 and K23 MH073091-01 (to RAM); National Health and Medical Research Council APP1173025 (to KLG).VA Career Development Award #1IK2CX002107 -US Department of Veterans Affairs CSR&D.ENIGMA was supported partly by NIH U54 EB020403 from the Big Data to Knowledge (BD2K) program, R56AG058854, R01MH116147, and P41 EB015922 (to PMT); NIMH R01MH106595 (to CMN).We thank Cohen Veterans Bioscience for ongoing support and building a collaborative scienti c environment.

Figures Figure 1
Figures