The genetic architecture of human cortical folding.

The folding of the human cerebral cortex is a highly genetically regulated process that allows for a much larger surface area to fit into the cranial vault and optimizes functional organization. Sulcal depth is a robust yet understudied measure of localized folding, previously associated with multiple neurodevelopmental disorders. Here, we report the first genome-wide association study of sulcal depth. Through the multivariate omnibus statistical test (MOSTest) applied to vertex-wise measures from 33,748 U.K. Biobank participants (mean age, 64.3 years; 52.0% female), we identified 856 genome-wide significant loci (P < 5 × 10−8). Comparisons with cortical thickness and surface area indicated that sulcal depth has higher locus yield, heritability, and effective sample size. There was a large amount of genetic overlap between these traits, with gene-based analyses indicating strong associations with neurodevelopmental processes. Our findings demonstrate sulcal depth is a promising neuroimaging phenotype that may enhance our understanding of cortical morphology.

INTRODUCTION

During early brain development, the cerebral cortical sheet folds into gyri and sulci in a highly regulated manner because of multiple intrinsic and extrinsic mechanical forces (–). This cortical folding not only allows for a much larger surface area to fit into the cranial vault but also reduces distance between neurons, leading to faster signal transmission (). Accordingly, measures of sulcal morphology are associated with cognitive performance (), and lack of cortical folding (lissencephaly) is accompanied by severe mental retardation (). Atypical folding can result from defects in neuronal proliferation, migration, and differentiation, and has been associated with major neurodevelopmental (–) and neurodegenerative disorders ().

Sulcal depth is a rather understudied measure of sulcal morphology, reflecting the convexity or concavity of any given point on the cortical surface. This measure is therefore highly suited for quantification of the primary folding pattern of the cortical sheet, and it has been shown to be insensitive to noise in the form of small wrinkles and to be relatively stable across individuals (). Sulcal depth is also known to have high test-retest reliability across the brain (, ). The few studies using this measure have provided evidence of regional differences between men and women () and have shown that sulci become more shallow with aging (, ), yet this remains to be fully characterized.

Several studies have indicated that there is a strong genetic component to sulcal depth, which is mostly prenatally determined (, , ). Sulci are more similar in monozygotic than in dizygotic twins (), and an estimated 56% of between-subject variance in average depth of the central sulcus is under genetic control (). Furthermore, Williams syndrome, caused by deletion of a section of chromosome 7, is associated with widespread reductions in sulcal depth (), which mediate its behavioral symptoms (, ). Yet, there has been no large-scale molecular genetics study of this measure.

Here, we provide the first genome-wide association study (GWAS) of sulcal depth, comparing its genetic architecture to the more commonly studied brain morphological measures of cortical thickness and surface area. Given that gene variants are likely to have distributed effects across magnetic resonance imaging (MRI) phenotypes, we targeted a multivariate analysis of a vertex-wise representation of the cortical surface, preventing the need for a multiple-comparison correction or data reduction strategies (). We therefore applied the multivariate omnibus statistical test (MOSTest) () to data from 1153 vertices, using a common template (fsaverage3) with the medial wall vertices excluded. Our primary sample consisted of 33,748 unrelated White British participants of the UK Biobank (UKB), with a mean age of 64.3 years (SD, 7.5 years); 52.0% were female. See figs. S1 and S2 for vertex-wise brain maps showing the distribution of the sulcal depth metric in this sample as well as its Pearson’s correlation with cortical thickness and surface area. We further carried out analyses in two additional samples consisting of (i) 5199 UKB participants (mean age, 62.8 years) that were excluded from the main analyses as they were not of White British descent, and (ii) 8072 participants of the Adolescent Brain Cognitive Development (ABCD) study (mean age, 9.9 years). These samples differed substantially from the discovery sample in terms of ethnicity and age, providing a strong test of generalizability of the reported associations. Following surface reconstruction, we preresidualized all vertices for age, sex, scanner site, a proxy of image quality (), and the first 20 genetic principal components to control for population stratification. After applying a rank-based inverse normal transformation, MOSTest was performed on the resulting residualized measures, yielding a multivariate association with each single-nucleotide polymorphism (SNP). We additionally repeated the main GWAS analyses while covarying for the mean across all vertices to remove global effects; the findings from these analyses were highly similar to the main analyses, as reported in the Supplementary Materials.

RESULTS

In the discovery sample, MOSTest revealed 856 independent loci, reaching the genome-wide significance threshold of α = 5 × 10−8 for sulcal depth (see Fig. 1A). In comparison, for surface area and thickness, we found 661 and 591 loci, respectively (see figs. S3 and S4). We found that the sets of discovered loci for sulcal depth had acceptable rates of generalization (see fig. S4 and table S2).

Fig. 1.

Locus discovery.

(A) Manhattan plot of the multivariate GWAS on sulcal depth, with the observed −log10(P value) of each SNP shown on the y axis. The x axis shows the relative genomic location, grouped by chromosome, and the red dashed line indicates the whole-genome significance threshold of 5 × 10−8. The y axis is clipped at −log10(P value) = 150. (B) Lateral view of the cortex, depicting the color-coded vertex-wise z values for the top hit rs4924345 at chromosome 15, across the discovery sample (top) and the two replication samples (middle and bottom). The left column corresponds to the left hemisphere, and the right column to the right hemisphere. (C) Power plot showing the relation between variance explained by genome-wide significant hits (y axis) and sample size (x axis). The number of hits discovered per modality and the percent explained genetic variance are indicated between brackets in the legend. (D) Bar plot of the mean SNP-based heritability (with 95% confidence interval) across vertices (on the y axis) per modality (x axis). In (C) and (D), sulcal depth is represented in green, surface area in orange, and cortical thickness in purple. ****P ≤ 0.0001.

The effects of discovered top variants on sulcal depth followed gyral and sulcal patterns, and the spatial pattern of effects was highly similar in the two additional samples, as shown in Fig. 1B for the most significant SNP at chromosome 15. Data files S1 to S3 contain information on all discovered loci per trait, including mapped genes, and list the significance of each lead SNP in each sample.

Next, using the MiXeR tool (, ), we fitted a Gaussian mixture model of the null and non-null effects to the three GWAS summary statistics, estimating the polygenicity and effect size variance (“discoverability”). The results are summarized in Fig. 1C, depicting the estimated proportion of genetic variance explained by discovered SNPs for each trait as a function of sample size. The horizontal shift of the curve across the different traits indicates that the effective sample size is the highest for sulcal depth and lowest for cortical thickness. Furthermore, the mean heritability of sulcal depth, calculated through linkage disequilibrium score regression (LDSC) () applied to each vertex, was significantly higher than for the two other traits (see Fig. 1D), i.e., the higher genetic signal in sulcal depth is also captured by univariate measures.

We additionally performed gene-based analyses through multimarker analysis of genomic annotation (MAGMA) (). We identified 2010 multiple comparison-corrected (α = 0.05/18,203) significant genes for sulcal depth (i.e., 11% of all 18,203 genes) versus 1486 for area and 1347 for thickness. The lists of all significant genes are provided in the Supplementary Materials.

Next, we analyzed the genetic overlap between the three traits at the locus level, gene level, and pathway level. At the locus level, we found that sulcal depth and surface area had 625 loci with overlapping start and end genomic positions (Dice coefficient of 0.82), 509 loci overlapped between sulcal depth and cortical thickness (Dice = 0.70), and 450 loci overlapped between surface area and thickness (Dice = 0.72). A total of 447 loci overlapped across all three traits (see Fig. 2A). We further found that a substantial number of the loci discovered for these cortical traits overlap with those previously identified for a set of brain disorders and cognitive performance (see table S3).

Fig. 2.

Genetic overlap.

(A) Venn diagram of the number of discovered loci overlapping between the three different traits. (B) Scatterplot of gene-based P values, with y axis indicating P values for sulcal depth and x axis for those for surface area, and the coloring indicating P values for cortical thickness. Note that −log10(P values) are clipped at 40. (C) Ten most significant gene pathways for sulcal depth, as listed on the y axis, with the −log10(P values) indicated on the x axis for each of the three traits. (D) Mean-normalized expression (y axis) of genes over time (x axis; log10 scale) per trait and over all available genes, as indicated by colors. Gray shading indicates 95% confidence bands.

The large genetic overlap between the traits was also evident at the gene level, as illustrated in Fig. 2B. The top gene STH, thought to play a role in phosphorylation of tau (), was highly significantly associated with all three traits. ROBO2, NAV2, and SEMA3A, key players in neuronal outgrowth guidance (–), were also associated with all three traits. The two histone genes, HIST1H4L and HIST1H2BL, both located in the large histone gene cluster at chromosome 6, were relatively specifically associated with sulcal depth; among its many roles, histone activity is central in regulating gene expression patterns that determine neuronal proliferation and differentiation processes that shape the cerebral cortex ().

Figure 2C shows results from gene set analyses through MAGMA (), listing the top 10 most significant gene ontology pathways for sulcal depth together with its P values for the two other traits. We found strong associations with neurogenesis and neuron differentiation pathways, overlapping between all three traits. Associations with neuronal tangential migration were shared by sulcal depth and surface area but much less by cortical thickness, in line with the role of tangential migration of neurons in determining cortical folding (). Notably, pathways related to chondrocyte differentiation and skeletal system development appeared more specific to sulcal depth, possibly pointing toward early life interactions between cortical folding and the shaping of the cranium (). Full results from the pathway analyses are provided in data files S7 to S9.

We further coupled the findings of our gene-based analyses to cortical gene expression patterns, derived from postmortem brain tissue of clinically unremarkable donors across the life span (). As shown in Fig. 2D, the probes tagging genes associated with the three traits showed a distinct profile over the life span, characterized by high prenatal expression and low postnatal expression.

DISCUSSION

Here, we reported the results from the first large-scale molecular genetics study of sulcal depth. With 856 loci discovered, explaining an estimated 32% of its genetic variance, this study has found the highest number of loci for any brain trait considered so far. We further provided evidence that our findings generalize to other populations and carry meaningful biological information.

The direct comparison with surface area and thickness indicated that sulcal depth is more heritable. This may reflect the evolutionary significance of cortical folding, the development of which enabled the advent of a larger brain and optimization of its functional organization (). A synthesis of the literature suggests that human-specific folding follows from an interplay between mechanical forces and cellular mechanisms that have come about over the course of evolution through mutations of genes primarily coupled to cell cycling and neurogenesis (). This is in accordance with the specific genes identified, known to play roles in the regulation of neuronal proliferation and migration (–). Our findings therefore suggest that the sulcal depth metric is closely aligned with these genetic processes that shape important brain morphological characteristics.

As indicated by the brain maps, genetic effects have opposing directions of effects on some neighboring points in the brain. This is in line with strong differences in the morphology and arrangement of neurons and fibers along cortical folds, varying widely from the gyral crown along the lateral wall down to the sulcal fundus (, ). Furthermore, the current findings indicate the presence of widespread genetic effects, illustrating this important characteristic of its genetic architecture; numerous variants are involved, each with a complex pattern of effects spread across the cortex. As shown by our estimates of generalization, these multivariate genetic effects on sulcal depth seem to be rather independent of ethnicity and age (), emphasizing the fundamental neurobiological relevance for brain morphology.

We further found large genetic overlap between all three morphological brain traits, extending our previous findings that surface area and thickness share the majority of their genetic determinants (, ). We found that the genetics of sulcal depth overlaps more with surface area than with cortical thickness, indicating a closer relation between the neurobiological mechanisms determining the degree of cortical folding and surface area. Patterns of tangential migration of neurons are thought to be important drivers of both folding of the cortical sheet and determination of its surface area (), while its thickness is influenced more by radial migration (). The results from the pathway analyses appear to confirm this distinction. However, the estimates of overlap, together with the observed phenotypic correlations, also indicate that sulcal depth and surface area do partly capture distinct genetic processes. Mammalian species also vary in their degree of folding independent of surface area (). Sulcal depth is therefore likely to provide additional information on the molecular genetic influences on brain morphology, capturing variation in the mechanical processes that determine the folding of the cortical sheet, to complement what is found through studies of surface area and cortical thickness. The identified genetic overlap with brain disorders and cognitive ability additionally indicates that the genetics of these brain measures carry clinically relevant information.

In addition to the reported locus overlap, the specific identified genetic variants, genes, and pathways further inform our understanding of cortical morphology and associated disorders. The most significant pathways were particularly relevant for early brain development, with neurogenesis and differentiation ranking highest. This fits very well with a large body of literature on the genetic regulation of the mechanical forces that drive cortical folding (). It is also in accordance with our findings that the sets of identified genes showed highest expression in fetal cortical tissue. Furthermore, cortical folding has been shown to take place almost entirely prenatally (), with sulcal patterns at birth being predictive of neurobehavioral outcomes (). Follow-ups on our work with neuroimaging data across the life span, including infants, are needed to replicate these findings and to further determine spatiotemporal patterns of genetic effects on sulcal depth. Given the reported associations of sulcal morphology with a range of neurodevelopmental and neurodegenerative disorders (–) as well as our exploratory findings of genetic overlap with brain disorders, it will also be of interest to investigate more thoroughly how sulcal depth genetics relates to the development of brain disorders over the life span.

We note as a limitation that we are currently lacking a strong mechanistic understanding of what the sulcal depth metric captures. Still, the evolutionary and ontogenetic importance of cortical folding appears to be reflected in the higher heritability compared to cortical thickness and surface area. Our exploration of multivariate genetic associations with sulcal depth attests to this, providing novel insight into the complexity of the human cortex, warranting further investigation into this measure to enhance our understanding of the brain and associated disorders.

MATERIALS AND METHODS

Participants

For the primary analyses, we made use of data from participants of the UKB population cohort, obtained from the data repository under accession number 27412. The composition, setup, and data gathering protocols of the UKB have been extensively described elsewhere (). UKB has received ethics approval from the National Health Service National Research Ethics Service (reference: 11/NW/0382) and obtained informed consent from its participants. For this study, we selected White British individuals, as determined by self-report and confirmed by genetic principal components analysis (UKB Data-Field 22006), which had undergone the neuroimaging protocol. We made use of T1 MRI scan data released up to March 2020, excluding 771 individuals with bad structural scan quality as indicated by an age- and sex-adjusted Euler number () more than 3 SDs lower than the scanner site mean. We further excluded one of each pair of related individuals, as determined through genome-wide complex trait analysis (GCTA), using a threshold of 0.0625 (n = 1138). Our sample size for this analysis was n = 33,748, with a mean age of 64.3 years (SD, 7.5). A total of 52.0% of the sample were female.

For the analyses of generalizability, we made use of the same UKB data and preprocessing steps but restricted our sample to those individuals who were not classified as White British. This left us with n = 5199 individuals, with a mean age of 62.8 years (SD, 7.7), of which 54.1% were female.

We additionally included data from children participating in the ABCD neurodevelopmental study, with complete genetic data and baseline T1 MRI scans from data release 3.0 [NIMH Data Archive (NDA) DOI:10.151.54/1519007] that passed the ABCD quality control procedures (n = 8072). These children had a mean age of 9.9 years (SD, 0.6), and 46.9% were female. All procedures were approved by a central Institutional Review Board (IRB) at the University of California, San Diego, and, in some cases, by individual site IRBs. Parents or guardians provided written informed consent, and children assented before participation.

Data preprocessing

UKB T1-weighted scans were collected from three scanning sites throughout the United Kingdom, all on identically configured Siemens Skyra 3T scanners, with a 32-channel receive head coils. The UKB core neuroimaging team has published extensive information on the applied scanning protocols and procedures, which we refer to for more details (). ABCD data were collected from 21 acquisition sites using Siemens Prisma, GE 750, and Philips 3T scanners to collect the T1-weighted scans. Scanning protocols were harmonized across sites. Full details of all imaging acquisition protocols used in ABCD are outlined by Casey et al. ().

All scans were stored locally at the secure computing cluster of the University of Oslo. We applied the standard “recon-all -all” processing pipeline of FreeSurfer v5.3, followed by extracting vertex-wise data for sulcal depth, surface area, and thickness, at ico3 (1284 vertices) and ico4 (5124) resolutions, without applying smoothing. We included both the left and right hemisphere measures and excluded noncortical vertices belonging to the medial wall.

Note that we have chosen sulcal depth as a metric of cortical folding, as it captures vertex-wise localized folding, providing the signed distance from the inflated surface.

We subsequently regressed out age, sex, scanner site, Euler number, and the first 20 genetic principal components from each vertex measure. Following this, we applied rank-based inverse normal transformation () to the residuals of each measure, leading to normally distributed measures as input for the GWAS.

We reran the MOSTest analyses as described above on the primary UKB sample, additionally regressing out the mean across all vertices for each of the three traits. The resulting number of loci is shown in table S1.

MOSTest procedure

The MOSTest software is freely available at https://github.com/precimed/mostest, and details about the procedure and its extensive validation have been described previously (). In brief, consider N variants and M (preresidualized) phenotypes. Let z be a z score from the univariate association test between ith variant and jth (residualized) phenotype, and z = (z, …, z) be the vector of z scores of the ith variant across M phenotypes. Let Z = {z} be the matrix of z scores, with variants in rows and phenotypes in columns. For each variant, consider a random permutation of its genotypes and let be the matrix of z scores from the univariate association testing between variants with permuted genotypes and phenotypes. A random permutation of genotypes is done once for each variant, and the resulting permuted genotype is tested for association with all phenotypes, therefore preserving correlation structure between phenotypes.

Let be the correlation matrix of , and is its singular valued decomposition (U and V are the orthogonal matrixes, and S is the diagonal matrix, with singular values of on the diagonal). Consider the regularized version of the correlation matrix , where S is obtained from S by keeping r largest singular values and replacing the remaining with rth largest. The MOSTest statistic for the ith variant (scalar) is then estimated as , where regularization parameter r is selected separately for cortical area and thickness to maximize the yield of genome-wide significant loci. In this study, we observed the largest yield for cortical surface area with r = 10; the optimal choice for cortical thickness was r = 20 and for sulcal depth r = 30. The distribution of the test statistics under null () is approximated from the observed distribution of the test statistics with permuted genotypes, using the empirical distribution in the 99.99th percentile and gamma distribution in the upper tail, where shape and scale parameters of gamma distribution are fitted to the observed data. The P value of the MOSTest test statistic for the ith variant is then obtained as .

Univariate GWAS procedure

We made use of the UKB v3 imputed data, which have undergone extensive quality control procedures as described by the UKB genetics team (). After converting the BGEN format to PLINK binary format, we additionally carried out standard quality check procedures, including filtering out individuals with more than 10% missingness, SNPs with more than 5% missingness, and SNPs failing the Hardy-Weinberg equilibrium test at P = 1 × 10−9. We further set a minor allele frequency threshold of 0.005, leaving 9,061,022 SNPs.

For ABCD, we took the genetic data that were part of the third data release, imputed through Trans-Omics for Precision Medicine (TOPMED), and mapped this back from genome build hg38 to hg19. We subsequently applied identical postimputation quality check procedures and filters for the UKB genetic data, leaving 13,131,314 SNPs.

We have previously calculated that the number of features provided by fsaverage3, 1153 vertices following exclusion of the medial wall, leads to the maximum number of loci identified through MOSTest, compared to other resolutions (). We therefore used the fsaverage3 resolution data for the input to MOSTest and for calculating heritability of the univariate vertex-wise data (described below). For visualization of the regionalization of the results, in the form of brain maps, we additionally carried out univariate GWAS on the 5124 vertices that make up fsaverage4, i.e., one level of resolution above fsaverage3. This was done only to improve the resolution of the visualizations. The univariate GWAS on each of the preresidualized and normalized measures were carried out using the standard additive model of linear association between genotype vector, g, and phenotype vector, y.

Heritability

We calculated the SNP-based heritability for each vertex at fsaverage3 resolution by applying LDSC to the univariate GWAS summary statistics with default settings ().

Locus definitions

Independent significant SNPs and genomic loci were identified from the MOSTest summary statistics in accordance with the psychiatric genomics consortium (PGC) locus definition, as also used in Functional Mapping and Annotation of GWAS (FUMA) SNP2GENE (). First, we select a subset of SNPs that pass the genome-wide significance threshold of 5 × 10−8 and use PLINK to perform a clumping procedure at linkage disequilibrium (LD) r2 = 0.6 to identify the list of independent significant SNPs. Second, we clump the list of independent significant SNPs at LD r2 = 0.1 threshold to identify lead SNPs. Third, we query the reference panel for all candidate SNPs in LD r2 of 0.1 or higher with any lead SNPs. Furthermore, for each lead SNP, its corresponding genomic loci is defined as a contiguous region of the lead SNPs’ chromosome, containing all candidate SNPs in r2 = 0.1 or higher LD with the lead SNP. Last, adjacent genomic loci are merged together if they are separated by less than 250 kb. Allele LD correlations are computed from the European (EUR) population of the 1000 Genomes phase 3 data. We additionally performed clumping according to the definition used by the Enhancing Neuroimaging Genetics through Meta-analysis Consortium to allow for comparison with previous imaging GWAS studies. According to this definition, loci were formed through PLINK using a P value threshold of 5 × 10−8 (--clump-p1) and LD cutoffs of 1 Mb (--clump-kb) and r2 < 0.2 (--clump-r2). Please see table S1 for the number of lead SNPs and loci according to both definitions. We made use of the FUMA online platform (https://fuma.ctglab.nl/) to map significant SNPs from the MOSTest analyses to genes.

Genetic overlap

Genetic overlap between the different MOSTest feature sets was operationalized as the number of significant loci that were physically overlapping between each pair of summary statistics. Loci from the pair are considered physically overlapping if their boundaries, their start and end genomic positions as determined through clumping, overlap. Contrary to metrics of global genetic correlation, this approach can be applied to MOSTest summary statistics and is insensitive to mixed directions of effects that would lower estimates of genetic correlation.

The Dice coefficient for each pair of traits was calculated as the number of overlapping loci divided by the sum of the total number of discovered loci for both traits.

MiXeR analysis

We applied a causal mixture model (, ) to estimate the percentage of variance explained by genome-wide significant SNPs as a function of sample size. For each SNP, i, MiXeR models its additive genetic effect of allele substitution, β, as a point normal mixture, , where π1 represents the proportion of non-null SNPs (“polygenicity”) and represents variance of effect sizes of non-null SNPs (discoverability). Then, for each SNP, j, MiXeR incorporates LD information and allele frequencies for 9,997,231 SNPs extracted from 1000 Genomes phase 3 data to estimate the expected probability distribution of the signed test statistic; , where N is sample size, H indicates heterozygosity of ith SNP, r indicates allelic correlation between ith and jth SNPs, and is the residual variance. Furthermore, the three parameters, , are fitted by direct maximization of the likelihood function. Fitting the univariate MiXeR model does not depend on the sign of z, allowing us to calculate ∣z∣ from MOSTest P values. Last, given the estimated parameters of the model, the power curve S(N) is then calculated from the posterior distribution p(δ∣z, N).

Gene set analyses

We carried out gene-based analyses using MAGMA v1.08 with default settings, which entail the application of an SNP-wide mean model and use of the 1000 Genomes phase 3 EUR reference panel to the three MOSTest summary statistics. Gene set analyses were carried out by applying MAGMA to the gene-level output, restricting the sets under investigation to those that are part of the Gene Ontology biological processes subset (n = 7522), as listed in the Molecular Signatures Database (c5.bp.v7.1).

Regarding the results from the gene-based analyses, in Fig. 2, we note that there is a horizontal line visible at –P = 5 × 10−10, caused by many genes having this exact P value. This is due to MAGMA switching to permutation when its numerical integration approach fails. MAGMA uses 1 × 10−9 permutations, so when the observed is more extreme than this, this is the resulting P value.

Gene expression analyses

We made use of gene expression data derived from brain tissue from 56 clinically unremarkable donors, ranging in age from 5 weeks post conception to 82 years (). We took the data as preprocessed in (), selecting for each gene the probe with the highest differential stability, n = 16,660. Given the relatively high homogeneity of expression patterns across cortical brain samples (), we subsequently averaged over 13 cortical regions, within donor, and normalized the expression values, within probe, across donors, to a range between 0 (lowest observed value) and 100 (highest observed value). As a check, we reran these analyses while restricting to samples from the individual cortical lobes. For each of the lobes, we observed expression patterns that were highly similar to those reported in the main analyses, with the gene sets showing high prenatal expression that dropped off at birth. Plotting of the mean expression over time per gene set was done with ggplot2 in R v4.0.3., with geom_smooth(method=”gam”) using default settings.