Genetic correlations and genome-wide associations of cortical structure in general population samples of 22,824 adults.

Cortical thickness, surface area and volumes vary with age and cognitive function, and in neurological and psychiatric diseases. Here we report heritability, genetic correlations and genome-wide associations of these cortical measures across the whole cortex, and in 34 anatomically predefined regions. Our discovery sample comprises 22,824 individuals from 20 cohorts within the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) consortium and the UK Biobank. We identify genetic heterogeneity between cortical measures and brain regions, and 160 genome-wide significant associations pointing to wnt/β-catenin, TGF-β and sonic hedgehog pathways. There is enrichment for genes involved in anthropometric traits, hindbrain development, vascular and neurodegenerative disease and psychiatric conditions. These data are a rich resource for studies of the biological mechanisms behind cortical development and aging.

Introduction

The cortex is the largest part of the human brain, associated with higher brain functions, such as perception, thought, and action. Brain cortical thickness (CTh), cortical surface area (CSA), and cortical volume (CV) are morphological markers of cortical structure obtained from magnetic resonance imaging (MRI). These measures change with age[1-3] and are linked to cognitive functioning[4,5]. The human cortex is also vulnerable to a wide range of disease or pathologies, ranging from developmental disorders and early onset psychiatric and neurological diseases to neurodegenerative conditions manifesting late in life. Abnormalities in global or regional CTh, CSA, and CV have been observed in neurological and psychiatric disorders, such as Alzheimer’s disease[6], Parkinson’s disease[7], multiple sclerosis[8], schizophrenia[9], bipolar disorder[9], depression[10], and autism[11]. The best method to study human cortical structure during life is using brain MRI. Hence, understanding the genetic determinants of the most robust MRI cortical markers in apparently normal adults could identify biological pathways relevant to brain development, aging, and various diseases. Neurons in the neocortex are organized in columns which run perpendicular to the surface of the cerebral cortex[12]; and, according to the radial unit hypothesis, CTh is determined by the number of cells within the columns and CSA is determined by the number of columns[13].

Thus, CTh and CSA reflect different mechanisms in cortical development[13,14] and are likely influenced by different genetic factors[15-18]. CV, which is the product of CTh and CSA, is determined by a combination of these two measures, but the relative contribution of CTh and CSA to CV may vary across brain regions. CTh, CSA, and CV are all strongly heritable traits[15-21] with estimated heritability of 0.69–0.81 for global CTh, and from 0.42 to 0.90 for global CSA[15,16,18]. Across different cortical regions, however, there is substantial regional variation in heritability of CTh, CSA, and CV[15-21].

Since CTh, CSA, and CV are differentially heritable and genetically heterogeneous, we explore the genetics of each of these imaging markers using genome-wide association analyses (GWAS) in large population-based samples. We study CTh, CSA, and CV in the whole cortex and in 34 cortical regions in 22,824 individuals from 21 discovery cohorts and replicate the strongest associations in 22,363 persons from the Enhancing Neuroimaging Genetics through Meta-analysis (ENIGMA) consortium. Our analyses reveal 160 genome-wide significant associations pointing to wnt/β-catenin, TGF-β, and sonic hedgehog pathways. We observe genetic heterogeneity between cortical measures and brain regions and find enrichment for genes involved in anthropometric traits, hindbrain development, vascular and neurodegenerative disease, and psychiatric conditions.

Results

Genome-wide association analysis

The analyses of global CTh, CSA, and CV included 22,163, 18,617, and 22,824 individuals, respectively. After correction for multiple testing (pDiscovery < 1.09 × 10−9), we identified no significant associations with global CTh. However, we identified 12 independent loci associated with global CSA (n = 6) and CV (n = 6). These are displayed in Supplementary Data 1 and Supplementary Figs. 1 and 2. Five of the 6 CSA loci were replicated in an external (ENIGMA consortium) sample[22]. The ENIGMA consortium only analyzed CSA and CTh.

GWAS of CTh, CSA, and CV in 34 cortical regions of interest (ROIs) identified 148 significant associations. There were 16 independent loci across 8 chromosomes determining CTh of 9 regions (Supplementary Data 2), 54 loci across 16 chromosomes associated with CSA of 21 regions (Supplementary Data 3), and 78 loci across 17 chromosomes determining CV of 23 cortical regions (Supplementary Data 4). We replicated 57 out of 64 regional CTh and CSA loci that were available in the ENIGMA consortium sample[22] using a conservative replication threshold of pReplication = 3.1 × 10−4, 0.05/160. Region-specific variants with the strongest association at each genomic locus are shown in Tables 1–3. Chromosomal ideograms showing genome-wide significant associations with global and regional cortical measures in the discovery stage are presented in Fig. 1.

Table 1

Genome-wide significant associations (pDiscovery < 1.09 × 10−9) of regional CTh.

Lobe	Region	Locus	Position	Lead SNP	Nearest gene	Annotation	N	pDiscovery	pReplication	ppooled
Temporal	Superior temporal	16q24.2	87225139	rs4843227	LOC101928708	Intergenic	21,887	2.79E−12	2.45E−05	2.31E−15
		17q21.31	44861003	rs199504	WNT3	Intronic	21,887	1.30E−10	1.17E−04	5.85E−13
	Middle temporal	14q23.1	59072144	rs10782438	KIAA0586	Intergenic	21,559	2.17E−13	2.76E−08	8.99E−21
	Inferior temporal	2q35	217332057	rs284532	SMARCAL1	Intronic	21,885	1.03E−09	2.64E−01	3.04E−07
	Banksts	14q23.1	59074878	rs160458	KIAA0586(	Intergenic	18,342	9.39E−10	2.42E−09	6.45E−18
Parietal	Superior parietal	16q24.2	87225101	rs9937293	LOC101928708	Intergenic	21,886	2.68E−14	1.64E−13	2.27E−27
		1q41	215141570	rs10494988	KCNK2	Intergenic	21,886	2.60E−12	3.66E−08	2.63E−19
	Postcentral	15q14	39633904	rs2033939	C15orf54	Intergenic	21,885	1.17E−73	5.18E−68	7.73E−136
Occipital	Lateral occipital	5q14.1	79933093	rs245100	DHFR	Intronic	21,886	2.68E−11	3.77E−06	1.16E−15
	Cuneus	14q23.1	59624317	rs4901904	DAAM1	Intergenic	21,885	4.02E−14	3.17E−10	2.88E−23
	Insula	16q12.1	51449978	rs7197215	SALL1	Intergenic	21,560	1.45E−13	2.00E−02	6.42E−12
		9q31.3	113679617	rs72748157	LPAR1	Intronic	21,560	1.46E−10	1.38E−04	5.16E−13

N number of individuals in meta-analysis, p two-sided p-value of discovery GWAS meta-analysis in CHARGE, p two-sided p-value of replication meta-analysis in ENIGMA, p two-sided p-value of pooled discovery and replication meta-analysis, p-values are not adjusted for multiple comparisons, banksts banks of the superior temporal sulcus.

in bold: significant replication—pReplication < 3.1 × 10−4 (= 0.05/Nl, Nl=160, total number of lead SNPs).

Table 3

Genome-wide significant associations (pDiscovery < 1.09 × 10−9) of global and regional CV.

Lobe	Region	Locus	Position	Lead SNP	Nearest gene	Annotation	N	pDiscovery
	Global	6q22.32	126792095	rs11759026	MIR588	Intergenic	22,410	6.31E−19
		17q21.31	44790203	rs169201	NSF	Intronic	22,784	2.11E−13
		17q21.32	43549608	rs149366495	PLEKHM1	Intronic	22,099	8.18E−13
		12q14.3	66358347	rs1042725	HMGA2	3’UTR	22,784	7.04E−11
		12q23.2	102921296	rs11111293	IGF1	Intergenic	22,784	5.45E−10
		6q22	109002042	rs4945816	FOXO3	3’UTR	22,784	8.93E−10
Frontal	Superior frontal	5q14.3	92186429	rs888814	NR2F1-AS1	Intergenic	22,692	3.29E−13
	Rostral middle frontal	15q14	39636227	rs17694988	C15orf54	Intergenic	22,793	3.15E−11
	Caudal middle frontal	2q12.1	105460333	rs745249	LINC01158	ncRNA_intronic	22,726	2.35E−11
		6q22.32	127068983	rs853974	RSPO3	Intergenic	22,351	4.82E−11
	Pars opercularis	5q23.3	128734008	rs12187568	ADAMTS19	Intergenic	20,753	4.27E−18
		15q14	39639898	rs4924345	C15orf54	Intergenic	22,758	1.97E−14
	Pars triangularis	3q24	147106319	rs2279829	ZIC4	UTR3	22,759	3.16E−23
		7q21.3	96196906	rs67055449	LOC100506136	Intergenic	22,759	4.03E−19
		15q14	39633904	rs2033939	C15orf54	Intergenic	22,759	8.49E−14
		7q21.3	96129071	rs62470042	C7orf76	Intronic	22,759	7.38E−13
		6q15	91942761	rs12660096	MAP3K7	Intergenic	22,759	4.74E−10
	Lateral orbitofrontal	14q22.2	54769839	rs6572946	CDKN3	Intergenic	22,801	2.29E−10
	Precentral	15q14	39634222	rs1080066	C15orf54	Intergenic	22,699	5.84E−125
		10q25.3	118648841	rs3781566	SHTN1	Intronic	22,699	4.68E−11
Temporal	Superior temporal	3q26.32	177296448	rs13084960	LINC00578	ncRNA_intronic	22,681	1.12E−11
	Banksts	14q23.1	59072226	rs186347	KIAA0586	Intergenic	22,727	1.15E−15
	Fusiform	14q23.1	59833172	rs1547199	DAAM1	Intronic	22,605	4.58E−10
		1p33	47980916	rs6658111	FOXD2	Intergenic	22,605	7.78E−10
	Transverse temporal	2q23.2	150012936	rs2046268	LYPD6B	Intronic	22,786	2.55E−12
	Parahippocampal	2q33.1	199809716	rs966744	SATB2	Intergenic	22,747	2.23E−10
Parietal	Superior parietal	15q14	39633904	rs2033939	C15orf54	Intergenic	22,723	4.28E−23
		16q24.2	87225139	rs4843227	LOC101928708	Intergenic	22,723	1.16E−13
		19p13.2	13109763	rs68175985	NFIX	Intronic	21,777	3.27E−11
		5q15	92866553	rs62369942	NR2F1-AS1	ncRNA_intronic	21,664	4.32E−10
	Inferior parietal	20q13.2	52448936	rs6097618	SUMO1P1	Intergenic	22,701	2.09E−17
		12q14.3	65797096	rs2336713	MSRB3	Intronic	22,701	2.47E−13
		3q13.11	104724634	rs971551	ALCAM	Intergenic	22,701	2.34E−10
	Supramarginal	15q14	39632013	rs71471500	THBS1	Intergenic	22,645	9.71E−28
		14q23.1	59627631	rs2164950	DAAM1	Intergenic	22,645	3.59E−20
		3q24	147106319	rs2279829	ZIC4	UTR3	22,645	5.36E−18
	Postcentral	15q14	39633904	rs2033939	THBS1	Intergenic	22,662	4.34E−133
		3q24	147106319	rs2279829	ZIC4	UTR3	22,662	2.54E−17
		9q21.13	76144318	rs67286026	ANXA1	Intergenic	22,662	5.03E−11
		2q36.3	226563259	rs16866701	NYAP2	Intergenic	22,545	5.69E−11
	Precuneus	14q23.1	59628609	rs74826997	DAAM1	Intergenic	22,803	4.85E−20
		3q28	190663557	rs35055419	OSTN	Intergenic	22,428	2.02E−10
		2p22.2	37818236	rs2215605	CDC42EP3	Intergenic	22,803	3.43E−10
		3q13.11	104713881	rs12495603	ALCAM	Intergenic	22,803	9.71E−10
Occipital	Lateral occipital	14q23.1	59627631	rs2164950	DAAM1	Intergenic	22,799	6.89E−16
	Lingual	14q23.1	59625997	rs73313052	DAAM1	Intergenic	22,805	1.06E−20
		6q22.32	127089401	rs2223739	RSPO3	Intergenic	22,805	1.75E−10
	Cuneus	14q23.1	59625997	rs73313052	DAAM1	Intergenic	22,799	4.59E−43
		11p15.3	12072213	rs11022131	DKK3	Intergenic	22,799	5.96E−12
		13q31.1	80192236	rs9545156	LINC01068	Intergenic	22,799	4.09E−10
	Pericalcarine	14q23.1	59628679	rs76341705	DAAM1	Intergenic	22,824	1.39E−29
		13q31.1	80191873	rs9545155	LINC01068	intergenic	22,824	2.25E−13
		11p14.1	30876113	rs273594	DCDC5	Intergenic	22,824	3.51E−13
		1p13.2	113208039	rs12046466	CAPZA1	Intronic	22,824	2.36E−12
		1p33	47980916	rs6658111	FOXD2	Intergenic	22,824	3.85E−11
		11q22.3	104012656	rs1681464	PDGFD	Intronic	22,824	7.51E−11
		6q22.32	127096181	rs9401907	RSPO3	Intergenic	22,824	2.11E−10
		7p21.1	18904400	rs12700001	HDAC9	Intronic	22,824	2.12E−10
		5q12.1	60315823	rs10939879	NDUFAF2	Intronic	22,824	2.92E−10
	Caudal anterior cingulate	5q14.3	82852578	rs309588	VCAN	Intronic	22,748	2.60E−10
	Insula	11q23.1	110949402	rs321403	C11orf53	Intergenic	22,543	9.58E−12
		8q24.12	120596023	rs10283100	ENPP2	Exonic	21,481	8.29E−11

Fig. 1

Chromosomal ideogram of genome-wide significant associations with measures of cortical structure.

Cortical surface areas, cortical volumes and cortical thickness. Each point represents the significantly associated variant, the colors correspond to the different cortical regions and the shape to different type of measument (pDiscovery < 1.09 × 10−9).

If we had used a more stringent threshold of pDiscovery < 4.76 × 10−10 = 5 × 10−8/105, correcting for all the 105 GWAS analyses performed, we would have identified 142 significant associations (Supplementary Data 1–4).

The strongest associations with CTh and CV were observed for rs2033939 at 15q14 (pDiscovery, CTh = 1.17 × 10−73 and pDiscovery, CV = 4.34 × 10−133) in the postcentral (primary somatosensory) cortex, and for CSA with rs1080066 at 15q14 (pDiscovery, CSA = 8.45 × 10−109) in the precentral (primary motor) cortex. Figure 2 shows the lowest p-value of each cortical region. The postcentral cortex was also the region with the largest number of independent associations, mainly at a locus on 15q14. The corresponding regional association plots are presented in Supplementary Fig. 3. Quantile-quantile plots of all meta-analyses are presented in Supplementary Figs. 4–7 and the corresponding genomic inflation factors (λGC), LD score regression (LDSR) intercepts, and ratios are shown in Supplementary Data 5. Although we observe inflated test statistics for some traits with λGC between 1.02 and 1.11, LDSR intercepts between 0.98 and 1.02 indicate that the inflation is mainly due to polygenicity. For traits with λGC > 1.05, the LDSR ratios range between 0.00 and 0.15 which means that a maximum of 15% of the inflation is due to other causes.

Fig. 2

Lowest discovery meta-analysis p-value of CSA, CTh, and CV in each cortical region.

a Lowest pDiscovery of CSA, b lowest pDiscovery of CTh, c lowest pDiscovery of CV.

Associations across cortical measures and with other traits

Supplementary Data 6 presents variants that are associated with the CSA or the CV across multiple regions. We observed 25 single nucleotide polymorphisms (SNPs) that determined both the CSA and CV of a given region, 4 SNPs that determined CTh and CV of the same region, but no SNPs that determined both the CTh and CSA of any given region (Supplementary Data 7). We also checked the overlap between our findings and two previous GWAS studies, including 8428[23] and 19,621[24] individuals from the UK Biobank, which among other phenotypes, investigate CTh, CSA, and CV (Supplementary Data 8). Regarding CTh, one variant, rs2033939 at 15q14, was associated with CTh of the postcentral gyrus in both studies. For CSA and CV, we found 11 associations at 15q14, 14q23.1 and 3q24, and 14 associations at 15q14, 14q23.1, 3q24, 8q24.1, 12q14.3, and 20q13.2, respectively, with the same cortical region as in our study. Out-of-sample polygenic risk score (PRS) analyses showed associations (pPRS < 4.76 × 10−3) with all investigated cortical measures in all cortical regions in 7800 UK Biobank individuals (Supplementary Data 9). For CTh, we observed the maximum phenotypic variance explained by the PRS (RPRS2) in the global cortex (RPRS2 = 0.015, pPRS = 1.05 × 10−26), and for CSA and CV in the pericalcarine cortex (RPRS2,CSA = 0.029, pPRS,CSA = 1.29 × 10−50; RPRS2,CV = 0.032, pPRS,CV = 5.30 × 10−56). When assessing genetic overlap with other traits, we observed that SNPs determining these cortical measures have been previously associated with anthropometric (height), neurologic (Parkinson’s disease, corticobasal degeneration, and Alzheimer’s disease), psychiatric (neuroticism and schizophrenia) and cognitive performance traits as well as with total intracranial volume (TIV) on brain MRI (Supplementary Data 10–12).

Gene identification

Positional mapping based on ANNOVAR showed that most of the lead SNPs were intergenic and intronic (Fig. 3). One variant, rs2279829, which was associated with both CSA and CV of the pars triangularis, postcentral and supramarginal cortices, is located in the 3′UTR of ZIC4 at 3q24. We also found an exonic variant, rs10283100, in gene ENPP2 at 8q24.12 associated with CV of the insula.

Fig. 3

Functional annotation categories for global and regional CTh, CSA, and CV.

Proportion of functional annotation categories for global and regional cortical thickness (blue), surface area (light green), and volume (yellow) assigned by ANNOVAR.

We used multiple strategies beyond positional annotation to identify specific genes implicated by the various GWAS associated SNPs. FUMA identified 232 genes whose expression was determined by these variants (eQTL) and these and other genes implicated by chromatin interaction mapping are shown in Supplementary Data 13–15. MAGMA gene-based association analyses revealed 70 significantly associated (p < 5.87 × 10−8) genes (Supplementary Data 16–18). For global CSA and CV, 7 of 9 genes associated with each measure overlapped, but there was no overlap with global CTh. For regional CSA and CV, we found 28 genes across 13 cortical regions that determined both measures in the same region. Figure 4 summarizes the results of GTEx eQTL, chromatin interaction, positional annotation, and gene-based mapping strategies for all regions. While there are overlapping genes identified using different approaches, only DAAM1 gene (Chr14q23.1) is identified by all types of gene mapping for CV of insula. eQTL associations of our independent lead SNPs in the Religious Orders Study Memory and Aging Project (ROSMAP) dorsolateral frontal cortex gene expression dataset are presented in Supplementary Data 19.

Fig. 4

Number of overlapping genes between gene mapping methods.

Number of overlapping genes between FUMA eQTL mapping, FUMA chromatin interaction mapping, ANNOVAR chromosome positional mapping, and MAGMA gene-based analysis for all cortical regions combined for cortical surface area (a), thickness (b) and volume (c).

Pathway analysis

MAGMA gene set analyses identified 7 pathways for CTh, 3 pathways for CSA and 9 pathways for CV (Supplementary Data 20). Among them are the gene ontology (GO) gene sets hindbrain morphogenesis (strongest association with thickness of middle temporal cortex), forebrain generation of neurons (with surface area of precentral cortex), and central nervous system neuron development (with volume of transverse temporal cortex). However, after Bonferroni correction only one significant pathway (p < 1.02 × 10−7) remained: regulation of catabolic process for CTh of the inferior temporal cortex. InnateDB pathway analyses of genes mapped to independent lead SNPs by FUMA showed a significant overlap between CTh and CSA genes and the Wnt signaling pathway (Supplementary Figs. 8 and 9) as well as a significant overlap between CV genes and the basal cell carcinoma pathway (Supplementary Fig. 10).

Heritability

Heritability estimates (h2) of global CTh were 0.64 (standard error (se) = 0.12; pSOLAR = 3 × 10−7) in the ASPS-Fam study and 0.45 (se = 0.08; pGCTA = 2.5 × 10−7) in the Rotterdam study (RS). For CSA, h2 was 0.84 (se = 0.12; pSOLAR = 2.63 × 10−11) in ASPS-Fam and 0.33 (se = 0.08, pGCTA = 1 × 10−4) in RS, and for CV, h2 was 0.80 (se = 0.11; pSOLAR = 1.10 × 10−9) in ASPS-Fam and 0.32 (se = 0.08; pGCTA = 1 × 10−4) in RS. There was a large range in heritability estimates of regional CTh, CSA, and CV (Supplementary Data 21).

Heritability based on common SNPs as estimated with LDSR was 0.25 (se = 0.03) for global CTh, 0.29 (se = 0.04) for global CSA and 0.30 (se = 0.03) for global CV. LDSR heritability estimates of regional CTh, CSA, and CV are presented in Supplementary Data 21 and Supplementary Fig. 11. For the regional analyses, the estimated heritability ranged from 0.05 to 0.18 for CTh, from 0.07 to 0.36 for CSA and from 0.06 to 0.32 for CV. Superior temporal cortex (h2CTh = 0.18, h2CSA = 0.30, h2CV = 0.26), precuneus (h2CTh = 0.16, h2CSA = 0.29, h2CV = 0.28) and pericalcarine (h2CTh = 0.15, h2CSA = 0.36, h2CV = 0.32) are among the most genetically determined regions.

The results of partitioned heritability analyses for global and regional CTh, CSA, and CV with functional annotation and additionally with cell-type-specific annotation are presented in Supplementary Data 22 and 23. For global CTh, we found enrichment for super-enhancers, introns and histone marks. Repressors and histone marks were enriched for global CSA, and introns, super-enhancers, and repressors for global CV. For regional CSA and CV the highest enrichment scores (>18) were observed for conserved regions.

Genetic correlation

We found high genetic correlation (rg) between global CSA, and global CV (rg = 0.81, pLDSR = 1.2 × 10−186) and between global CTh and global CV (rg = 0.46, pLDSR = 1.4 × 10−14), but not between global CTh and global CSA (rg = −0.02, pLDSR = 0.82). Whereas the genetic correlation between CSA and CV was strong (rg > 0.7) in most of the regions (Supplementary Data 24 and Supplementary Fig. 12), it was generally weak between CSA and CTh with rg < 0.3, and ranged from 0.09 to 0.69 between CTh and CV. The postcentral and lingual cortices were the two regions with the highest genetic correlations between both CTh and CV, as well as CTh and CSA.

Genetic correlation across the various brain regions for CTh (Supplementary Fig. 13, Supplementary Data 25), CSA (Supplementary Fig. 14, Supplementary Data 26), and CV (Supplementary Fig. 15, Supplementary Data 27) showed a greater number of correlated regions for CTh and greater inter-regional variation for CSA and CV. Supplementary Data 28–30 and Supplementary Figs. 16–18 show genome-wide genetic correlations between the cortical measures and anthropometric, neurological and psychiatric, and cerebral structural traits.

Discussion

In our genome-wide association study of up to 22,824 individuals for MRI determined cortical measures of global and regional thickness, surface area, and volume, we identified 160 genome-wide significant associations across 19 chromosomes. Heritability was generally higher for cortical surface area and volume than for thickness, suggesting a greater susceptibility of cortical thickness to environmental influences. We observed strong genetic correlations between surface area and volume, but weak genetic correlation between surface area and thickness. We identified the largest number of novel genetic associations with cortical volumes, perhaps due to our larger sample size for this phenotype, which was assessed in all 21 discovery samples.

It is beyond the scope of our study to discuss each of the 160 associations identified. A large number of the corresponding genes are involved in pathways that regulate morphogenesis of neurons, neuronal cell differentiation, and cell growth, as well as cell migration and organogenesis during embryonic development. At a molecular level, the wnt/β-catenin, TGF-β, and sonic hedgehog pathways are strongly implicated. Gene-set-enrichment analyses revealed biological processes related to brain morphology and neuronal development.

Broad patterns emerged showing that genes determining cortical structure are also often implicated in development of the cerebellum and brainstem (KIAA0586, ZIC4, ENPP2) as well as the neural tube (one carbon metabolism genes DHFR and MSRBB3, the latter also associated with hippocampal volumes[25]). These genes determine development of not only neurons but also astroglia (THBS1) and microglia (SALL1). They determine susceptibility or resistance to a range of insults: inflammatory, vascular (THBS1, ANXA1, ARRDC3-AS1[26]) and neurodegenerative (C15orf53, ZIC4, ANXA1), and have been associated with pediatric and adult psychiatric conditions (THBS1).

There is a wealth of information in the supplementary tables that can be mined for a better understanding of brain development, connectivity, function and pathology. We highlight this potential by discussing in additional detail, the possible significance of 6 illustrative loci, 5 of which, at 15q14, 14q23.1, 6q22.32, 17q21.31, and 3q24, associate with multiple brain regions at low p-values, while the locus at 8q24.12 identifies a plausible exonic variant.

The Chr15q14 locus was associated with cortical thickness, surface area, and volumes in the postcentral gyrus as well as with surface area or volume across six other regions in the frontal and parietal lobes. Lead SNPs at this locus were either intergenic between C15orf53 and C15orf54, or intergenic between C15orf54 and THBS1 (Thrombospondin-1). C15orf53 has been associated with an autosomal recessive form of spastic paraplegia showing intellectual disability and thinning of the corpus callosum (hereditary spastic paraparesis 11, or Nakamura Osame syndrome). Variants of THBS1 were reported to be related to autism[27] and schizophrenia[28]. The protein product of THBS1 is involved in astrocyte induced synaptogenesis[29], and regulates chain migration of interneuron precursors migrating in the postnatal radial migration stream to the olfactory bulb[30]. Moreover, THBS1 is an activator of TGFβ signaling, and an inhibitor of pro-angiogenic nitric oxide signaling, which plays a role in several cancers and immune-inflammatory conditions.

Variants at Chr14q23.1 were associated with cortical surface area and volume of all regions in the occipital lobe, as well as with thickness, surface area, and volume of the middle temporal cortex, banks of the superior temporal sulcus, fusiform, supramarginal and precuneus regions, areas associated with discrimination and recognition of language or visual form. These variants are either intergenic between KIAA0586, the product of which is a conserved centrosomal protein essential for ciliogenesis, sonic hedgehog signaling and intracellular organization, and DACT1, the product of which is a target for SIRT1 and acts on the wnt/β-catenin pathway. KIAA0586 has been associated with Joubert syndrome, another condition associated with abnormal cerebellar development. Other variants are intergenic between DACT1 and DAAM1 or intronic in DAAM1. DAAM1 has been associated with occipital lobe volume in a previous GWAS[31].

Locus 6q22.32 contains various SNPs associated with cortical surface area and volume globally, and also within some frontal, temporal and occipital regions. The SNPs are intergenic between RSPO3 and CENPW. RSPO3 and CENPW have been previously associated with intracranial[32,33] and occipital lobe volumes[31]. RSPO3 is an activator of the canonical Wnt signaling pathway and a regulator of angiogenesis.

Chr17q21.31 variants were associated with global cortical surface area and volume and with regions in temporal lobe. These variants are intronic in the genes PLEKHM1, CRHR1, NSF, and WNT3. In previous GWAS analyses, these genes have been associated with general cognitive function[34] and neuroticism[35]. CRHR1, NSF, and WNT3 were additionally associated with Parkinson’s disease[36] and intracranial volume[32,33,37]. The NSF gene also plays a role in Neuronal Intranuclear Inclusion Disease[38] and CRHR1 is involved in anxiety and depressive disorders[39]. This chromosomal region also contains the MAPT gene, which plays a role in Alzheimer’s disease, Parkinson’s disease, and frontotemporal dementia[40,41].

The protein product of the gene ZIC4 is a C2H2 zinc finger transcription factor that has an intraneuronal, non-synaptic expression and auto-antibodies to this protein have been associated with subacute sensory neuronopathy, limbic encephalitis, and seizures in patients with breast, small cell lung or ovarian cancers. ZIC4 null mice have abnormal development of the visual pathway[42] and heterozygous deletion of the gene has also been associated with a congenital cerebellar (Dandy-Walker) malformation[43], thus implicating it widely in brain development as well as in neurodegeneration. C2H2ZF transcription factors are the most widely expressed transcription factors in eukaryotes and show associations with responses to abiotic (environmental) stressors. Another transcription factor, FOXC1, also associated with Dandy-Walker syndrome has been previously shown to be associated with risk of all types of ischemic stroke and with stroke severity. Thus, ZIC4 might be a biological target worth pursuing to ameliorate neurodegenerative disorders.

We found an exonic SNP within the gene ENPP2 (Autotaxin) at 8q24.12 to be associated with insular cortical volume. This gene is differentially expressed in the frontal cortex of Alzheimer patients[44] and in mouse models of Alzheimer disease, such as the senescence-accelerated mouse prone 8 strain (SAMP8) mouse. Autotaxin is a dual-function ectoenzyme, which is the primary source of the signaling lipid, lysophosphatidic acid. Besides Alzheimer disease, changes in autotaxin/lysophosphatidic acid signaling have also been shown in diverse brain-related conditions, such as intractable pain, pruritus, glioblastoma, multiple sclerosis, and schizophrenia. In the SAMP8 mouse, improvements in cognition noted after administration of LW-AFC, a putative Alzheimer remedy derived from the traditional Chinese medicinal prescription ‘Liuwei Dihuang’ decoction, are correlated with restored expression of four genes in the hippocampus, one of which is ENPP2.

Among the other genetic regions identified, many have been linked to neurological and psychiatric disorders, cognitive functioning, cortical development, and cerebral structure (detailed listing in Supplementary Data 31).

Heritability estimates are, as expected, generally higher in the family-based Austrian Stroke Prevention-Family study (ASPS-Fam) than in the Rotterdam Study (RS) for CTh (average h2ASPS-Fam = 0.52; h2RS = 0.26), CSA (0.62 and 0.30) and CV (0.57 and 0.23). This discrepancy is explained by the different heritability estimation methods: pedigree-based heritability in ASPS-Fam versus heritability based on common SNPs that are in LD with causal variants[45] in RS.

Average heritability over regions is also higher for surface area and volume, than for thickness. The observed greater heritability of CSA compared to CTh is consistent with the previously articulated hypothesis, albeit based on much smaller numbers, that CSA is developmentally determined to a greater extent with smaller subsequent decline after young adulthood, whereas CTh changes over the lifespan as aging, neurodegeneration and vascular injuries accrue[1,3]. It is also interesting that brain regions more susceptible to early amyloid deposition (e.g., superior temporal cortex and precuneus) have a higher heritability.

We found no or weak genetic correlation between CTh and CSA, globally and regionally, and no common lead SNPs, which indicates that these two morphological measures are genetically independent, a finding consistent with prior reports[15,16]. In contrast, we found strong genetic correlation between CSA and CV and identified common lead SNPs for CSA and CV globally, and in 12 cortical regions. Similar findings have been reported in a previous publication[16]. The genetic correlation between CTh and CV ranged between 0.09 and 0.77, implying a common genetic background in some regions (such as the primary sensory postcentral and lingual cortices), but not in others. For CTh, we observed genetic correlations between multiple regions within each of the lobes, whereas for CSA and CV, we found genetic correlations mainly between different regions of the occipital lobe. Chen et al.[46] have also reported strong genetic correlation for CSA within the occipital lobe. There were also a few genetic correlations observed for regions from different lobes, suggesting similarities in cortical development transcended traditional lobar boundaries.

A limitation of our study is the heterogeneity of the MR phenotypes between cohorts due to different scanners, field strengths, MR protocols and MRI analysis software. This heterogeneity as well as the different age ranges in the participating cohorts may have caused different effects over the cohorts. We nevertheless combined the data of the individual cohorts to maximize the sample size as it has been done in previous CHARGE GWAS analyses[31-33]. To account for the heterogeneity we used a sample-size weighted meta-analysis that does not provide overall effect estimates. This method has lower power to detect associations compared to inverse-variance weighted meta-analysis and we therefore may have found less associations.Our inability to replicate 8 of the 76 genome-wide significant findings for CTh and CSA could be caused by false-positive results but may also be explained by insufficient power due to a too small sample size. Moreover, our sample comprises of mainly European ancestry, limiting the generalizability to other ethnicities. Strengths of our study are the population-based design, the large age range of our sample (20–100 years), use of three cortical measures as phenotypes of cortical morphometry, and the replication of our CTh and CSA findings in a large and independent cohort. In conclusion, we identified patterns of heritability and genetic associations with various global and regional cortical measures, as well as overlap of MRI cortical measures with genetic traits and diseases that provide new insights into cortical development, morphology, and possible mechanisms of disease susceptibility.

Methods

Study population

The sample of this study consist of up to 22,824 participants from 20 population-based cohort studies collaborating in the Cohorts of Heart and Aging Research in Genomic Epidemiology (CHARGE) consortium and the UK Biobank (UKBB). All the individuals were stroke- and dementia free, aged between 20 and 100 years, and of European ancestry, except for ARIC AA with African ancestry. Supplementary Data 32 provides population characteristics of each cohort and the Supplementary Methods provide a short description of each study. Each study secured approval from institutional review boards or equivalent organizations, and all participants provided written informed consent. Our results were replicated using summary GWAS findings of 22,635 individuals from the ENIGMA consortium.

Genotyping and imputation

Genotyping was conducted using various commercially available genotyping arrays across the study cohorts. Prior to imputation, extensive quality control was performed in each cohort. Genotype data were imputed to the 1000 Genomes reference panel (mainly phase 1, version 3) using validated software. Details on genotyping, quality control and imputation can be found in Supplementary Data 33.

Phenotype definition

This study investigated CTh, CSA, and CV globally in the whole cortex and in 34 cortical regions. Global and regional CTh was defined as the mean thickness of the left and the right hemisphere in millimeter (mm). Global CSA was defined as the total surface area of the left and the right hemisphere in mm2, while regional CSA was defined as the mean surface area of the left and the right hemisphere in mm2. Global and regional CV was defined as the mean volume of the left and the right hemisphere in mm3. The 34 cortical regions are listed in the Supplementary Methods. High resolution brain magnetic resonance imaging (MRI) data was obtained in each cohort using a range of MRI scanners, field strengths and protocols. CTh, CSA, and CV were generated using the Freesurfer software package[47] in all cohorts except for FHSucd, where an in-house segmentation method was used. MRI protocols of each cohort can be found in Supplementary Data 35 and descriptive statistics of CTh, CSA, and CV can be found in Supplementary Data 36–38.

Based on a predefined analysis plan, each study fitted linear regression models to determine the association between global and regional CTh, CSA, and CV and allele dosages of SNPs. Additive genetic effects were assumed and the models were adjusted for sex, age, age[2], and if needed for study site and for principal components to correct for population stratification. Cohorts including related individuals calculated linear mixed models to account for family structure. Details on association software and covariates for each cohort are shown in Supplementary Data 33. Models investigating regional CTh, CSA, and CV were additionally adjusted for global CTh, global CSA and global CV, respectively. Quality control of the summary statistics shared by each cohort was performed using EasyQC[48]. Genetic variants with a minor allele frequency (MAF) <0.05, low imputation quality (R2 < 0.4), and which were available in less than 10,000 individuals were removed from the analyses. Details on quality control are provided in the Supplementary Methods.

We then used METAL[49] to perform meta-analyses using the z-scores method, based on p-values, sample size, and direction of effect, with genomic control correction. To estimate the number of independent tests for the p-value threshold correction, we used a non-parametric permutation testing procedure[50-53] in the combined Rotterdam Study cohort (N = 4442) and UK Biobank (N = 8213). First, we generated a random independent variable, to insure that there is no true relationship between brain measurements and this variable. Second, we ran linear regression analyses between this variable and all brain measurements one-by-one in each of the cohorts separately (104 regressions in total per cohort). Third, we saved the minimum p-value obtained from those 104 regressions. Then, as suggested in literature[54], we repeated this procedure 10.000 times. Therefore, at the end we had 10.000 minimum p-values per cohort. The minimum p-value distribution follows a Beta distribution Beta(m,n), where m = 1 and n is the degree of freedom, which represents the number of independent tests in case of permutation testing. Using python statistical library we fitted the Beta function with the saved minimum p-values, and found n for Rotterdam Study and UK Biobank identically equal to 46. Based on the permutation test results, the genome-wide significance threshold was set a priori at 1.09 × 10−9 (=5 × 10−8/46). We used the clumping function in PLINK[55] (linkage disequilibrium (LD) threshold: 0.2, distance: 300 kb) to identify the most significant SNP in each LD block. We used LDSR to calculate genomic inflation factors (λGC), LDSR intercepts and LDSR ratios for each meta-analysis. The LDSR intercept was estimated to differentiate between inflation due to a polygenic signal and inflation due to population stratification[56]. The LDSR ratio represents the amount of inflation that is due to other causes than polygenicity such as population stratification or cryptic relatedness.

For replication of our genome-wide significant CTh and CSA associations, we used GWAS meta-analysis results from the ENIGMA consortium[22] for all SNPs that were associated at a p-value <5 × 10−8 and performed a pooled meta-analysis. The p-value threshold for replication was set to 3.1 × 10−4 (=0.05/160: nominal significance threshold divided by total number of lead SNPs). CV was not available in the ENIGMA results. PRS analysis was performed for 7800 out of sample subjects (not included in the current GWAS) from UK Biobank cohort using the PRSice-2 software[57] with standard settings. The significance threshold for the association between the PRS and the phenotype was set to 4.76 × 10−3 (=0.05/105: nominal significance divided by number GWAS phenotypes). The NHGRI-EBI Catalog of published GWAS[58] was searched for previous SNP-trait associations at a p-value of 5 × 10−8 of lead SNPs. Regional association plots were generated with LocusZoom[59], and the chromosomal ideogram with PHENOGRAM (http://visualization.ritchielab.org/phenograms/plot).

Annotation of genome-wide significant variants was performed using the ANNOVAR software package[60] and the FUMA web application[61]. FUMA eQTL mapping uses information from three data repositories (GTEx, Blood eQTL browser, and BIOS QTL browser) and maps SNPs to genes based on a significant eQTL association. We used a false discovery rate threshold (FDR) of 0.05 divided by number of tests (46) to define significant eQTL associations. Gene-based analyses, to combine the effects of SNPs assigned to a gene, and gene set analyses, to find out if genes assigned to significant SNPs were involved in biological pathways, were performed using MAGMA[62] as implemented in FUMA. The significance threshold was set to 5.87 × 10−8 (=0.05/18522*46: FDR threshold divided by number of genes and independent tests) for gene-based analyses and to 1.02 × 10−7 (=0.05/10651: FDR threshold divided by the number of gene sets) for the gene set analyses. Additionally, FUMA was used to investigate a significant chromatin interaction between a genomic region in a risk locus and promoter regions of genes (250 bp upstream and 500 bp downstream of a TSS). We used an FDR of 1 × 10−6 to define significant interactions.

We investigated cis (<1 Mb) and trans (>1 MB or on a different chromosome) expression quantitative trait loci (eQTL) for genome-wide significant SNPs in 724 post-mortem brains from ROSMAP[63,64] stored in the AMP-AD database. The samples were collected from the gray matter of the dorsolateral prefrontal cortex. The significance threshold was set to 0.001 (=0.05/46: FDR threshold divided by the number of independent tests). For additional pathway analyses of genes that were mapped to independent lead SNPs by FUMA, we searched the InnateDB database[65]. The STRING database[66] was used for visualizing protein–protein interactions. Only those protein subnetworks with five or more nodes are shown.

Additive genetic heritability (h2) of CTh, CSA, and CV was estimated in two studies: the Austrian Stroke Prevention Family Study (ASPS-Fam; n = 365) and the Rotterdam Study (RS, n = 4472). In the population-based family study ASPS-Fam, the ratio of the genotypic variance to the phenotypic variance was calculated using variance components models in SOLAR[67]. In case of non-normalty, phenotype data were inverse-normal transformed. In RS, SNP-based heritability was computed with GCTA[68]. These heritability analyses were adjusted for age and sex.

Heritability and partitioned heritability based on GWAS summary statistics was calculated from GWAS summary statistics using LDSR) implemented in the LDSC tool (https://github.com/bulik/ldsc). Partitioned heritability analysis splits genome-wide SNP heritability into 53 functional annotation classes (e.g., coding, 3′UTR, promoter, transcription factor binding sites, conserved regions etc.) and additionally to 10 cell-type specific classes (e.g., central nervous system, cardiovascular, liver, skeletal muscle, etc.) as defined by Finucane et al.[69] to estimate their contributions to heritability. The significance threshold was set to 2.05 × 10−5 (=0.05/53*46: nominal significance divided by number of functional annotation classes and number of independent tests) for heritability partitioned on functional annotation classes and 2.05 < 10−6 (=0.05/53*10*46: nominal significance divided by number of functional annotation classes, number of cell types and number of independent tests) for heritability partitioned on annotation classes and cell types.

LDSR genetic correlation[70] between CTh, CSA, and CV was estimated globally and within each cortical region. The significance threshold was set to 7.35 × 10−4 (nominal threshold (0.05) divided by number of regions (34) and by number of correlations (CSA and CV, CSA and CTh). Genetic correlation was also estimated between all 34 cortical regions for CTh, CSA, and CV, with the significance threshold set to 8.91 × 10−5 (nominal threshold (0.05) divided by number of regions (34) times the number of regions −1 (33) divided by 2 (half of the matrix). Additionally, the amount of genetic correlation was quantified between CTh, CSA, and CV and physical traits (height, body mass index), neurological and psychiatric diseases (e.g., Alzheimer’s disease, Parkinson’s disease), cognitive traits and MRI volumes (p-value threshold (0.05/46/number of GWAS traits). As recommended by the LDSC tool developers, only HapMap3 variants were included in these analyses, as these tend to be well-imputed across cohorts.

Table 2

Genome-wide significant associations (pDiscovery < 1.09 × 10−9) of global and regional CSA.

Lobe	Region	Locus	Position	Lead SNP	Nearest gene	Annotation	N	pDiscovery	pReplication	ppooled
	Global	17q21.31	44787313	rs538628	NSF	Intronic	18,617	1.78E−23	4.45E−22	1.00E−43
		6q22.32	126792095	rs11759026	MIR588	Intergenic	18,617	5.21E−22	1.45E−14	3.50E−34
		6q22.33	127204623	rs9375477	RSPO3	Intergenic	18,617	4.86E−13	1.60E−08	1.23E−19
		6q21	109000316	rs9398173	FOXO3	Intronic	18,617	6.84E−10	2.96E−03	2.05E−10
Frontal	Superior frontal	5q14.3	92187932	rs17669337	NR2F1-AS1	Intergenic	18,272	1.40E−11	2.05E−06	8.07E−16
	Caudal middle frontal	6q22.32	126876580	rs9388500	RSPO3	Intergenic	17,891	2.35E−11	NA	NA
	Pars opercularis	5q23.3	128734008	rs12187568	ADAMTS19	Intergenic	16,632	1.19E−16	NA	NA
	Pars triangularis	3q24	147106319	rs2279829	ZIC4	UTR3	18,265	6.32E−20	1.94E−27	1.20E−45
		7q21.3	96175094	rs10458281	LOC100506136	Intergenic	18,265	1.15E−17	2.42E−11	1.20E−26
	Precentral	15q14	39634222	rs1080066	C15orf54	Intergenic	18,267	8.45E−109	2.53E−95	1.00E−200
		6q15	92002569	rs9345124	MAP3K7	Intergenic	18,267	5.50E−11	2.73E−14	9.91E−24
Temporal	Superior temporal	2p16.3	48274592	rs386645843	FBXO11	Intergenic	18,269	9.51E−12	8.42E−07	1.71E−16
		4q26	119249835	rs55699931	PRSS12	Intronic	18,269	2.08E−11	2.72E−02	6.96E−10
		2q23.2	150022681	rs13008194	LYPD6B	Intronic	18,269	5.94E−11	2.54E−07	1.92E−16
	Middle temporal	6q22.32	126964510	rs4273712	RSPO3	Intergenic	18,269	6.93E−10	1.07E−04	1.99E−12
	Banksts	14q23.1	59072226	rs186347	KIAA0586	Intergenic	18,265	4.11E−10	1.83E−09	4.93E−18
	Fusiform	17q21.31	44822662	rs199535	NSF	Intronic	18,269	1.01E−13	1.14E−06	8.13E−18
	Transverse temporal	2q23.2	150012936	rs2046268	LYPD6B	Intronic	18,264	9.09E−10	3.21E−10	1.78E−18
Parietal	Superior parietal	15q14	39632013	rs71471500	C15orf54	Intergenic	18,270	3.85E−24	5.55E−19	5.88E−41
		19p13.2	13109763	rs68175985	NFIX	Intronic	17,324	8.84E−11	2.68E−17	2.90E−26
	Inferior parietal	20q13.2	52448936	rs6097618	SUMO1P1	Intergenic	18,267	1.78E−16	NA	NA
		12q14.3	65797096	rs2336713	MSRB3	Intronic	18,267	1.24E−12	2.99E−12	2.85E−23
		2p25.2	4563477	rs669952	LINC01249	Intergenic	18,267	4.47E−10	1.37E−08	4.73E−17
	Supramarginal	15q14	39633904	rs2033939	C15orf54	Intergenic	18,272	9.07E−27	1.61E−28	1.59E−53
		14q23.1	59627631	rs2164950	DAAM1	Intergenic	18,272	1.25E−13	3.79E−14	3.46E−26
		3q24	147106319	rs2279829	ZIC4	UTR3	18,272	7.38E−12	4.24E−16	2.29E−26
	Postcentral	15q14	39634222	rs1080066	C15orf54	Intergenic	18,265	5.65E−47	2.44E−36	1.87E−80
		3q24	147106319	rs2279829	ZIC4	UTR3	18,265	1.90E−21	1.69E−26	2.92E−46
		9q21.13	76144318	rs67286026	ANXA1	Intergenic	18,265	3.58E−12	8.04E−06	7.82E−16
	Precuneus	14q23.1	59628609	rs74826997	DAAM1	Intergenic	18,270	2.40E−24	4.41E−18	4.59E−40
		6q23.3	138866268	rs9376354	NHSL1	Intronic	18,270	7.80E−13	4.12E−08	7.28E−19
		3q26	190666643	rs1159211	SNAR-I	Intergenic	18,270	4.49E−10	2.04E−05	1.59E−13
Occipital	Lateral occipital	14q23.1	59627631	rs2164950	DAAM1	Intergenic	18,269	3.04E−26	2.92E−15	2.25E−38
	Lingual	14q23.1	59628679	rs76341705	DAAM1	Intergenic	18,270	1.57E−20	8.67E−13	9.96E−31
	Cuneus	14q23.1	59625997	rs73313052	DAAM1	Intergenic	18,267	1.90E−32	3.19E−15	2.96E−43
		13q31.1	80191873	rs9545155	LINC01038	Intergenic	18,267	5.15E−10	2.98E−05	3.91E−13
	Pericalcarine	14q23.1	59628679	rs76341705	DAAM1	Intergenic	18,267	4.67E−24	2.56E−19	3.35E−41
		5q12.1	60117723	rs6893642	ELOVL7	Intronic	18,267	1.40E−13	1.68E−08	6.29E−20
		3q13.11	104724787	rs971550	ALCAM	Intergenic	18,267	2.18E−10	1.31E−06	4.49E−15
		6q22.33	127185801	rs9375476	RSPO3	Intergenic	18,267	2.20E−10	2.24E−08	4.32E−17
		1p13.2	113239478	rs2999158	MOV10	Intronic	18,267	6.46E−10	8.39E−10	3.49E−18
		13q31.1	80191873	rs9545155	LINC01068	Intergenic	18,267	7.51E−10	7.53E−09	4.05E−17
	Posterior cingulate	5q12.3	66104105	rs17214309	MAST4	Intronic	18,268	7.84E−11	1.52E−05	4.04E−14
	Insula	10q25.3	118704077	rs1905544	SHTN1	Intronic	17,599	4.06E−12	3.65E−03	1.28E−11

N number of individuals in meta-analysis, p two-sided p-value of discovery GWAS meta-analysis in CHARGE, p two-sided p-value of replication meta-analysis in ENIGMA, p two-sided p-value of pooled discovery and replication meta-analysis, p-values are not adjusted for multiple comparisons, banksts banks of the superior temporal sulcus.

NA, SNP or region not available in the ENIGMA sample.

In bold: significant replication—pReplication < 3.1  × 10−4 (= 0.05/Nl, Nl = 160, total number of lead SNPs).