The genetic architecture of structural left–right asymmetry of the human brain

1The human brain is characterized by various population-level asymmetries on its left–right axis, including an overall ‘torque’ whereby the left hemisphere extends posteriorly and ventrally relative to the right, a left–right difference in frontal–occipital gradients of cortical thickness, and hemispheric differences of morphology around the Sylvian fissure. Many brain functions are also lateralized, including hand motor control and language which show left-hemisphere dominance in roughly 85% of people–. Altered brain or behavioural asymmetries have been reported in various cognitive and psychiatric disorders,–, which suggests that population-typical asymmetries are linked to optimal human brain function.

2Behavioural and anatomical brain asymmetries are already apparent in utero,–, which indicates an early genetic-developmental programme of brain left–right axis formation,. Studies of visceral organ development (the heart, stomach, liver, etc.) have revealed that the generation of population-level asymmetry requires at least three important steps in the early embryo,: (1) the breaking of bilateral symmetry to create a left–right axis in a consistent orientation relative to the anterior–posterior and dorsal–ventral axes, (2) the triggering of different patterns of gene expression on the left and right sides of early embryonic structures and (3) the translation of asymmetric gene expression into lateralized morphology and organ placement.

3In principle, establishing an embryonic left–right axis requires chirality at some level, that is, key biomolecules or cellular structures that exist in only one of two possible mirror forms. Life on Earth is based on L-form amino acids rather than the mirror D-form, and this chirality carries through to the macrostructure and movement of primary cilia,, which help to create the left–right axis of the visceral organs in embryos. However, hemispheric dominances for language and hand motor control do not typically reverse in people with situs inversus of the viscera, that is, reversal of the visceral organs on the left–right axis, when caused by mutations in genes that encode primary ciliary components–. This observation indicates that there are distinct and possibly organ-intrinsic mechanisms at play in brain development, but such mechanisms remain unidentified. Therefore the genetic origins of human brain asymmetry remain unknown.

4Post mortem studies that contrasted gene expression between the left and right sides of the embryonic central nervous system have identified molecular pathways that may be involved, including that the two sides may transition through developmental stages slightly out of synchrony with each other–. However, these studies were necessarily based on data from just handfuls of samples, because of limited availability arising from appropriate ethical, legal and practical concerns. An alternative approach to identify genetic influences on brain asymmetry is to relate genomic variation in large population datasets to variation in adult brain asymmetry. Only three loci have previously been reported at a genome-wide significant level to affect variation in adult human brain asymmetries, in studies targeted at single features of temporal lobe anatomy, and they did not yield broader insights into biological pathways,. Here we made use of the unprecedented sample size available through the 2020 release of the UK Biobank magnetic resonance imaging (MRI) data, in combination with genome-wide genotype data, to perform the first multivariate, brain-wide genetic analysis of human brain anatomical asymmetry. The results from this analysis also allowed us to test whether genetic polymorphisms that are associated with variation in brain asymmetry are also associated with neurodevelopmental disorders or other behavioural and psychological traits, using publicly available genome-wide association scan (GWAS) summary statistics for these traits.

5For each of 32,256 participants with post-quality-control MRI and genetic data (Methods), and each of 73 bilaterally paired regional measures of brain structure, we calculated hemispheric asymmetry indexes (AI) as (left − right)/((left + right)/2) (for 33 cortical surface area AIs, 33 cortical thickness AIs and 7 subcortical volume AIs; Supplementary Table 1). The measures were derived from cortical parcellation and subcortical segmentation of structural brain images (Methods). All but one of the regional mean AIs were significantly different from zero, indicating population-level asymmetries (Bonferroni-corrected P < 0.05, Supplementary Fig. 1 and Supplementary Table 2), consistent with previous reports,. For example, some language-related regions showed greater average left than right surface areas, including superior temporal and supramarginal cortex, and pars opercularis.

6GCTA software was used to estimate the single-nucleotide polymorphism (SNP)-based heritability (h2) for each AI, that is, the extent to which variance in each AI was linked to common genetic variation over the entire genome (Methods). Forty-two AIs showed significant SNP-based heritabilities (false discovery rate (FDR)-corrected P < 0.05), that is, 28 of the surface area AIs, 8 cortical thickness AIs and 6 subcortical volume AIs (Fig. 1a and Supplementary Table 3), ranging from 2.2% for the AI of entorhinal cortical thickness to 9.4% for the AI of superior temporal surface area. The overall pattern was consistent with previous, twin-based heritability analyses,.

7SNP-based genetic correlation analysis (again using GCTA software) indicated overlapping genetic contributions to some of the AIs (Fig. 1b, Supplementary Fig. 2 and Supplementary Tables 4–10). Within a small number of cortical regions (Supplementary Table 10), surface area and thickness AIs had negative genetic correlations, which indicates that variants can have antagonistic effects on surface and thickness asymmetries of these regions.

8We performed a multivariate (mv)GWAS for 9,803,522 SNPs, using meta-canonical correlation analysis as implemented in MetaPhat, with the 42 AIs that had significant SNP-based heritability. This analysis tested each SNP separately for its simultaneous associations with all 42 AIs (Methods). A multivariate approach had the dual advantages of achieving data reduction and increasing statistical power compared with running 42 separate univariate GWAS. FUMA software was used to clump mvGWAS results on the basis of linkage disequilibrium (LD), and identify lead SNPs (maximally associated SNPs) at each associated locus. There were 21 distinct genomic loci at the 5 × 10−8 significance level associated with different aspects of brain asymmetry (Fig. 2, Table 1 and Supplementary Fig. 3), represented by 27 independent lead SNPs (with pairwise LD r2 < 0.1, Table 1).

9For each lead SNP, phenotype decomposition identified the ‘central’ AIs that contributed to its multivariate association (Supplementary Table 11). Most central AIs associated with the 27 lead SNPs were distributed in core regions of the language system (for example, lateral temporal, pars opercularis and supramarginal) or limbic system (for example, cingulate, orbitofrontal and mesial temporal cortex; Fig. 3). For example, the most significant SNP, with multivariate r = 0.094, P = 4.75 × 10−38 (rs41298373 on 10p14) had five central AIs: the minor allele was associated with a leftward shift of surface area asymmetry for two lateral temporal regions, a rightward shift of surface area asymmetry for two medial temporal regions and a leftward shift of cortical thickness asymmetry in the inferior temporal gyrus (Supplementary Table 11). All lead SNPs were associated with at least one cortical regional surface area AI as one of their central traits. Thirteen lead SNPs were associated with cortical thickness AIs, and five lead SNPs were associated with subcortical volume AIs. Two lead SNPs, rs35853889 and rs6658111, were associated with the AIs of both surface area and cortical thickness within the same region, that is, parahippcampal area and thickness AIs were associated with both of these SNPs, and rostral anterior cingulate area and thickness AIs were associated with rs35853889. In addition, five lead SNPs affecting cortical surface area AIs were also associated with AIs of subcortical volumes (Table 1 and Fig. 3), and the locus on 17q21 was associated with asymmetries of cortical surface area, thickness and subcortical volume (Table 1). The univariate associations of lead variants separately with left and right hemispheric measures are presented in Supplementary Table 11.

10FUMA software applied three strategies to annotate candidate SNPs to genes at significantly associated loci (Methods): physical position, expression quantitative trait locus (eQTL) information and chromatin interactions (Supplementary Table 12 and Supplementary Figs. 4 and 5). Here we summarize and cite notable information on each of the lead SNPs:

11Ten of the loci had annotations involving cytoskeleton-related genes: On 1p33, rs6658111 is close to AL356458.1, a pseudogene of MTMR14 (myotubularin related protein 14). On 2p23.3, rs62130503 is intronic to MAPRE3 (microtubule associated protein RP/EB family member 3a), and rs12617392 is a brain eQTL of MAPRE3. Located between these two SNPs on 2p23.3 is also AGBL5, which is a post-translational modifier of tubulin,. On 2q34, rs368536282 is close to MAP2 (microtubule associated protein 2), a well-known dendrite-specific marker of neurons, previously implicated in left-handedness by GWAS analysis–. On 6p21.33, rs7781 is in the 3′ untranslated region (UTR) of TUBB (tubulin beta class I). On 12q13.12, rs10783306 is close to the alpha tubulin gene TUBA1B. This variant is also in high LD with a handedness-associated variant, rs11168884 (r2 = 0.89). On 14q23.1, two lead variants for two independent genomic loci, rs160459 and rs201816193, showed evidence for cross-locus chromatin interaction via the promoters of nearby genes in foetal cortex (Supplementary Fig. 4). The former is near to DACT1, a locus which has been reported to associate with superior temporal sulcus depth, while the latter is close to DAAM1, which modulates the reorganization of the actin cytoskeleton and the stabilization of microtubules,. Two lead variants on 16q24.3, rs72813426 and rs111398992, are in introns of SPIRE2 and the tubulin gene TUBB3, respectively, both of which are key proteins in cytoskeleton organization,. rs111398992 is also in high LD with a handedness-associated variant, rs4550447 (r2 = 0.94). On 17q21.31 there were five independent lead SNPs: rs35908989 is intronic to MAPT which encodes microtubule-associated protein tau, and rs55938136, rs35853889 and rs568039055 are brain eQTLs,, of MAPT, while rs80103986 is in high LD (r2 = 0.91) with handedness-associated variant rs55974014. On 19p13.3, rs11672092 is intronic to the tubulin gene TUBB4A, and in high LD with rs66479618 (r2 = 0.88), another handedness-associated SNP.

12The 11 other loci did not have obvious microtubule-related annotations, but most had annotations related to brain phenotypes or development: On 3q24, rs2279829 is a cortical eQTL of ZIC4, which is involved in visual and auditory pathway development. rs9307052 on 4q22.1 is in high LD (r2 = 0.99) with the handedness-associated variant rs28658282. On 5q15, rs869219775 is close to NR2F1, which is involved in neural activity during cortical patterning. On 6q22.31-q22.32, rs9385385 is close to NCOA7, a nuclear receptor co-activator with its most abundant expression in the brain. On 7p14.3, rs6947352 is intronic to BBS9, which causes Bardet–Biedl syndrome when mutated, involving retinopathy and intellectual disability,. On 9q22.33, rs911934 is located in a region having a chromatin interaction with TRIM14 in adult cortex (Supplementary Fig. 4), a gene which may activate Wnt/β-catenin signalling and affects mesodermal versus ectodermal differentiation of embryonic stem cells. On 10p14, rs41298373 is a predicted deleterious missense coding variant in ITIH5, which was previously reported to affect planum temporale volumetric asymmetry. Inter-alpha-trypsin inhibitor proteins are involved in extracellular matrix stabilization. On 20p12.1, rs6135555 is in a region having a chromatin interaction with the FLRT3 promoter in neural progenitor cells (Supplementary Fig. 4), a gene which regulates axon guidance and excitatory synapse development. On 21q22.3, rs7283026 is intronic to COL18A1, involved in neural tube closure and mutated in Knobloch syndrome, which can include skull abnormalities. On 22q13.31, rs9615351 is an exonic variant of a gene involved in planar cell polarity, CELSR1. On Xp22.33, rs12400461 is close to pseudogene ASS1P4 and upstream of MXRA5; the latter encodes a matrix remodelling-associated protein and is implicated in autism.

13To further link asymmetry-associated SNPs to genes, the MAGMA software was used to perform genome-wide gene-based association analysis based on the results from mvGWAS. There were 57 significant genes at Bonferroni-corrected P < 0.05 (Supplementary Fig. 6 and Supplementary Table 13). Five of these genes were previously associated with handedness: MAP2, FAM13A, TUBA1B, TUBB3 and CRHR1. Forty-three of the 57 genes have been reported to associate with educational attainment and 15 with intelligence (Supplementary Table 14). For the proteins encoded by the 57 genes, there were 80 known or putative pairwise interactions in the STRING database, compared with 8 interactions expected for a random set of this size from the whole proteome (P < 1 × 10−16). This observation supports the validity of our mvGWAS association findings, as random noise would not lead to such functional clustering. Microtubule-related genes (for example, MAP2, MAPT, SPIRE2 and TUBA1A) linked different clusters together in the largest protein interaction network (Fig. 4a).

14We also used the genome-wide, gene-based P values for functional enrichment analysis using MAGMA, in relation to 7,343 Gene Ontology ‘biological process’ sets defined within the MSigDB database. The gene set ‘regulation_of_microtubule_binding’ (P = 3.73 × 10−6) showed significant enrichment (adjusted P < 0.05, Bonferroni correction, Supplementary Table 15). Significant enrichment within various microtubule-related sets, such as ‘microtubule_cytoskeleton_organization’ (P = 2.19 × 10−7) and ‘microtubule_based_process’ (P = 2.36 × 10−6), was also found when using the list of single closest genes (Table 1) to the 27 lead SNPs (Supplementary Table 16). Enrichment in microtubule-related sets was not reported in a recent GWAS of bilaterally averaged cortical surface area and thickness measures in 51,665 individuals, suggesting a particular involvement in hemispheric asymmetry rather than bilateral measures. We observed no statistically significant relation of our gene-based association P values with differential expression across cell types (Methods).

15Testing our genome-wide, gene-based P values with respect to human gene expression data from the BrainSpan database, from either 29 age groups or 11 defined developmental stages, we found relatively higher mRNA expression of brain-asymmetry-associated genes during early-prenatal (P = 4.27 × 10−3) and early–mid-prenatal (P = 9.37 × 10−4) stages, from 9 (P = 1.84 × 10−3) to 24 (P = 7.36 × 10−3) post-conceptional weeks (FDR-corrected P values <0.05) (Fig. 4b,c and Supplementary Table 17). This is consistent with the fact that various anatomical asymmetries of the brain are already visible in utero,, and supports the existence of an early developmental mechanism for establishing the brain’s left–right axis,,.

16We next used iSECA software to perform genetic overlap analysis with our mvGWAS results in relation to GWAS summary statistics from neurodevelopmental disorders, behavioural and psychological traits which have been reported to associate phenotypically with aspects of structural and/or functional brain asymmetry: attention deficit hyperactivity disorder,–, autism spectrum disorder,–, educational attainment,,, handedness,,, intelligence,– and schizophrenia,–. There was evidence for genetic overlap between brain asymmetries and autism (P = 0.005), educational attainment (P = 0.001) and schizophrenia (P = 0.002) which remained significant at Bonferroni-corrected P < 0.05 (Fig. 5, Supplementary Figs. 7 and 8 and Supplementary Table 18). In other words, SNPs that showed lower (more significant) association P values in our mvGWAS for brain asymmetry showed a statistically significant tendency to also show lower P values in previous, large-scale GWAS of autism, educational attainment and schizophrenia. Although we did not observe genetic overlap of brain asymmetry with handedness at a genome-wide level, we did note individual loci in common between these traits (above). In addition, we found no overlap between our mvGWAS results and those from a previous GWAS of intracranial volume in 32,438 participants (Supplementary Table 18, Supplementary Figs. 7 and 8), which again indicates that the genetic architecture of brain asymmetry is largely distinct from brain size.

17To achieve a reasonable level of genetic homogeneity in our mvGWAS, we had excluded any individuals not annotated as having ‘white British ancestry’ (through a combination of self-report and clustering based on principal components that capture major axes of ancestral diversity from the genotype data (Methods)). The UK Biobank includes additional participants who self-identify as being ‘white’, but who did not self-identify as British, or did not cluster genetically with the bulk of the ‘white British’ ancestry participants (Methods). After applying the same quality control criteria to these additional participants as in our discovery mvGWAS (except with respect to ancestry), and imposing the extra criterion that relatedness kinship coefficients should not be greater than 0.0442 with any participants from the discovery mvGWAS, data were available from 3,600 participants for an independent replication set. We tested each of the 27 lead SNPs from the mvGWAS in the replication set, using the same approach as the mvGWAS analysis, except that 40 genetic principal components were used as covariates to control for the greater degree of ancestral diversity in the replication set.

18Ten of the 27 independent lead SNPs from the discovery mvGWAS showed association P values < 0.05 on multivariate testing in the replication set (Supplementary Table 19). The combined P value of the remaining 17 lead SNPs was P = 3.3 × 10−4 in the replication set (calculated by Stouffer’s method), which we confirmed by permutation with respect to 10,000 repeat random samplings of 17 SNPs from the whole genome in the replication set (permutation-based P = 4 × 10−4). This indicates that the limited sample size of the replication set, compared with the discovery set, did not provide adequate power to replicate at the level of some individual SNPs, but that in combination there was evidence for replication. Moreover, among the 17 SNPs that showed P > 0.05 on multivariate testing in the replication set, some showed association P < 0.05 on univariate testing of the specific central traits identified for those SNPs in the discovery mvGWAS (Supplementary Table 19). It is also worth noting that 4 of these 17 SNPs (or SNPs in high linkage disequilibrium with them) have been reported to associate with left-handedness at a genome-wide significant level (see above for details), which is an additional form of validation with respect to a phenotype related to brain asymmetry. As also mentioned above, the high degree of functional clustering of genes identified through gene-based association testing is another form of support for the association results in the mvGWAS.

19Despite the importance of asymmetry as an organizing feature of the human brain, the early developmental processes which establish its left–right axis are unknown. In this study, we carried out multivariate GWAS analysis which identified 21 genetic loci associated with different aspects of adult brain anatomical asymmetry. Functional annotation implicated genes particularly involved in microtubule organization and prenatal brain development. Our study therefore sheds new light on the molecular genetic foundations of human brain asymmetry. In addition, at the genome-wide level, SNPs associated with brain asymmetry overlapped with those associated with educational attainment, autism and schizophrenia, while five specific loci that are associated with both brain asymmetry and handedness were identified.

20Previous studies in invertebrates and frog embryos have shown that the cytoskeleton plays a role in determining cellular chirality and asymmetrical patterning of other organs–. Cellular chirality refers to directional biases in cellular morphology, position, rotation or migration, which arise because of the inherent chirality of intracellular macromolecules such as those composing the cytoskeleton, and also manifest in terms of the intracellular distributions of organelles. For example, during early cell divisions in Xenopus (frog) embryos, the cytoskeleton has been reported to mediate asymmetric intracellular protein localization, as protein transport molecules move along cytoskeletal tracks within cells. Thus the cytoskeleton can provide a directionally consistent, organ-intrinsic bias during embryonic development that acts as a determinant of morphological asymmetry,,, arising from fundamental aspects of molecular and cellular biology.

21As our study associated cytoskeletal and embryonically expressed genes with variation in adult human brain asymmetry, it is possible that these genes are involved in the establishment of left–right asymmetry during early brain development, through a mechanism involving cellular chirality. As mentioned in Main, at least some aspects of human brain asymmetry appear to be uncoupled from the developmental pathway that leads to left–right organization of the visceral organs, which involves cilia and the nodal pathway. A cytoskeletal-mediated mechanism of brain asymmetry may therefore be organ intrinsic,,, that is, distinct from other pathways that establish broader aspects of body asymmetry.

22In this study, we identified genetic loci that are associated with 42 heritable aspects of brain asymmetry through a multivariate, brain-wide approach. A multivariate approach can boost statistical power while achieving data reduction, compared with separate univariate analyses of individual brain traits. A single set of genome-wide association results, pertaining simultaneously to multiple aspects of brain asymmetry, was then taken forward into functional annotation and downstream analyses, such as testing for genetic overlaps with other traits. The multivariate approach therefore helped to detect and interpret key aspects of the genetic architecture of brain asymmetry, without the noise inherent in repeat univariate testing. An important challenge remained to identify the particular brain traits that drove the multivariate associations at each locus, which was achieved in MetaPhat by decomposing associations into sets of ‘central’ traits on the basis of the Bayesian information criterion and canonical correlation P values.

23A consequence of the multivariate approach is that it does not yield univariate association effect sizes, and therefore mvGWAS results cannot be used for standard genetic correlation analyses, such as is performed with LD score regression. Therefore, we used iSECA to explore the genetic overlap of brain asymmetry with neurodevelopmental disorders, behavioural and psychological traits. This analysis was based on SNP-wise P values of association over the whole genome, that is, a genetic overlap was found when the SNPs tending to show low P values in our mvGWAS for brain asymmetry also tended to show low P values in publicly available GWAS summary statistics for a given disorder or trait. We found significant genetic overlaps of brain asymmetry with autism, schizophrenia and educational attainment, which suggests that genes affecting brain asymmetry also influence these traits. This is in line with literature that has shown phenotypic associations between altered brain asymmetry and these traits (see Main), and indicates that such phenotypic associations are contributed to an extent by genetic factors. As we found that brain asymmetry-related genes tend to be especially highly expressed in the embryonic brain, it seems likely that the genetic overlap of brain asymmetry and disorders reflects a genetic susceptibility to alterations of early neurodevelopment away from the typical trajectory. However, brain asymmetry continues to develop throughout the lifespan,, and the UK Biobank consists of middle- to older-age adults, so that our mvGWAS may have also identified genetic factors that affect brain asymmetrical changes later in life.

24Further research, for example, using Mendelian randomization, will be needed to understand whether brain asymmetries mediate gene–disorder associations in a causal sense, or whether altered brain asymmetry and disorder susceptibility are two distinct consequences arising from a partly overlapping genetic basis. It will also be important to map, on a brain-regional basis, which aspects of asymmetry show the strongest genetic overlaps with disorder susceptibility. Larger imaging–genetic datasets may be needed to support causal mediation and mapping studies with respect to disorders, as the present genetic association analysis was based on brain-wide asymmetry (albeit in a multivariate context).

25Many of the brain asymmetries were strong and directional at the population level, but their heritabilities were generally low, ranging up to 10%. This suggests that developmental mechanisms for brain asymmetry are tightly constrained and largely genetically invariant in the population, and that environmental factors and/or developmental randomness are responsible for most variability,–. A cytoskeleton-based origin of brain asymmetry would fit this scenario, as the cytoskeleton is essential for various fundamental functions in cellular biology, beyond axis formation,. Previous, twin- and family-based analyses, have reported heritabilities up to roughly 25% for some of the same asymmetry measures that we analysed in the present study, with an overall similar regional pattern. Higher heritabilities were found particularly for regions that are important in the language system (for example, superior temporal cortex) and limbic system (for example, medial temporal and cingulate cortex). Twin-based heritability is often measured to be higher than SNP-based heritability, which may be expected because SNPs are just one class of genetic variation, and also because twin studies can overestimate heritability when certain assumptions are not fully met. As twin studies have not indicated effects of shared environment on brain aymmetries,, early developmental randomness is likely to cause most variation.

26We did not correct for handedness as a covariate in our genetic analyses, as it is generally not advisable to correct for covariates which are themselves partly heritable. This is because biased genetic effects can be measured with respect to the target trait (in this case, brain asymmetry). Handedness cannot therefore be treated safely as a confound variable when analysing brain asymmetry. We were interested in any genetic effects associated with brain asymmetry, regardless of whether they might also be shared with other traits such as handedness. Having identified genetic variants associated with brain asymmetry, we then queried post hoc whether they have been reported as significant in previous GWAS of handedness in over 1 million people. We did not observe a significant genetic overlap between structural brain asymmetry and handedness at the genome-wide level, which again may be due to the relatively low SNP-based heritabilities of these traits, in combination with limited statistical power in the present sample size for this kind of analysis. However, five individual SNPs associated with both brain asymmetry and handedness were identified, which suggests that a significant genome-wide overlap might be detected when using a larger dataset in the future.

27The UK Biobank currently includes by far the largest single imaging–genetic dataset available to the scientific community. A limitation of the present study is the lack of a large, age-matched replication sample with comparable homogeneity of ancestry to the discovery mvGWAS analysis. We included an independent replication sample of 3,600 individuals from the UK Biobank, with greater diversity of ancestry than the 32,256 individuals of the primary mvGWAS. Association was replicated for 10 of the 27 independent lead SNPs from the mvGWAS, and for the remaining 17, their combined P value in the replication set was 3.3 × 10−4, indicating that a larger dataset would have supported replication at the level of more individual SNPs. Functional clustering of the closest genes to the 27 independent lead SNPs, according to microtubule-related biology and protein–protein interactions, also supported the validity of the mvGWAS findings, as this was unlikely to occur by chance. In addition, four of the SNPs that did not replicate individually at P < 0.05 (or SNPs in high linkage disequilibrium with them; see Results) were previously associated at a genome-wide significant level with left-handedness, which is an additional form of validation with respect to a behavioural asymmetry trait.

28Another limitation of the present study is that functional enrichment analysis with respect to post-mortem, developmental gene expression was only possible from roughly 8 weeks post conception and onwards. The Brainspan database does not contain sufficient numbers of samples from earlier embryonic stages than this. It may be that the left–right brain axis is established extremely early in development, for example during formation of the neural tube, which begins in the third week post conception. As noted in Main, possibilities for research using human embryos from this developmental stage are highly restricted. Gene expression studies using animal models may therefore be necessary to understand the establishment of the mammalian brain’s left–right axis.

29In conclusion, our findings motivate genetic-developmental studies of left–right differentiation of the embryonic mammalian brain, focused on a possible cytoskeletal-mediated mechanism of axis formation. Our study also suggests that disruption of this mechanism may contribute to susceptibility to cognitive and psychiatric disorders, in line with asymmetry being an important aspect of healthy brain organization for many functions.

30This study was conducted under UK Biobank application 16066, with C.F. as principal investigator. The UK Biobank is a general adult population cohort. The UK Biobank received ethical approval from the National Research Ethics Service Committee North West-Haydock (reference 11/NW/0382), and all of their procedures were performed in accordance with the World Medical Association guidelines. Informed consent was obtained for all participants. We used the brain imaging data released in February 2020, and data availability and processing (described below) resulted in a final sample of 32,256 participants of white British ancestry for the primary GWAS analysis, together with the structural MRI data and genotype data from the same participants. The age range of these participants was 45 to 81 years (mean 63.77 years), 15,288 were male and 16,968 were female. An independent replication dataset of 3,600 individuals was also drawn from the UK Biobank, who self-identified as white, but not British, or did not cluster genetically with the bulk of the ‘white British’ ancestry participants (see below). The age range of these participants was 45 to 80 years (mean 62.89 years), 1,574 were male and 2,026 were female.

31We downloaded imputed SNP genotype data from the UK Biobank data portal (bgen files; imputed data v3-release March 2018). We first excluded subjects with a mismatch of their self-reported and genetically inferred sex, with putative sex chromosome aneuploidies, or who were outliers on the basis of heterozygosity (principle-component-corrected heterozygosity >0.19) and genotype missingness (missing rate >0.05) as calculated by Bycroft et al.. The primary analyses were restricted to participants with ‘white British ancestry’, which was defined by Bycroft et al. (‘in.white.British.ancestry.subset’), using a combination of self-report and cluster analysis on the basis of the first six principal components which capture genetic ancestry from the genome-wide genotype data. We randomly excluded one from each pair of remaining individuals who had a kinship coefficient >0.0442, as defined by Bycroft et al.. Next, QCTOOL (v.2.0.6) and PLINK were used to perform genotype quality control: excluding SNPs with minor allele frequency <1%, Hardy–Weinberg equilibrium test P value <1 × 10−7 and imputation INFO score <0.7 (a measure of genotype imputation confidence). We also excluded multi-allelic SNPs because most of the downstream software (below) could not handle them. This resulted in 9,803,522 bi-allelic variants.

32The same process was applied to derive the independent replication dataset of 3,600 individuals, except that these participants did not self-identify as British (although they did identify as ‘white’), or they did not fall within the bounds of the ‘white British ancestry’ cluster as defined by Bycroft et al.. We also imposed the extra requirement that participants in the replication set should not show a kinship coefficient >0.0442 with any individual in the primary discovery dataset.

33Brain anatomical measures of regional cortical surface area, cortical thickness and subcortical volumes were derived from the structural scans (Siemens Skyra 3-T MRI with 32-channel radiofrequency receive head coil) released by the UK Biobank Imaging Study (for the full protocol, see http://biobank.ndph.ox.ac.uk/showcase/refer.cgi?id=2367). Briefly, in vivo whole-brain T1-weighted MRI scans were used to perform cortical parcellation into 34 regions per hemisphere with the Desikan–Killiany atlas, and 7 subcortical structural segmentations. Surface area was measured at the grey–white matter boundary, and thickness was measured as the average distance in a region between the white matter and pial surfaces. Details of the image quality control and processing are described elsewhere. Given that the data for the temporal pole were reported as unreliable, we only used 33 surface area, 33 cortical thickness and 7 subcortical volume measures in each hemisphere (Supplementary Table 1). Per measure, we removed data points greater than six standard deviations from the mean. Then, we calculated the AI for each matching pair of left and right measures, in each participant, as (left − right)/((left + right)/2). Given this definition, a positive AI reflects leftward asymmetry (greater left than right). The AI is a widely used measure in brain asymmetry studies,,. The denominator ensures that the index does not simply scale with brain size, that is, the left–right difference is adjusted for the bilateral measure. For each AI, one-sample t testing was used to examine whether the population mean AI was significantly different from zero, with Bonferroni correction at 0.05 for multiple testing. Subsequently, the distributions of AIs were normalized by rank-based inverse normalization to minimize statistical artifacts. The normalized AIs were used as input for subsequent analysis.

34The Desikan–Killiany atlas was derived from manual segmentations of sets of reference brain images. The labelling system incorporates hemisphere-specific information on sulcal and gyral geometry with spatial information regarding the locations of brain structures, and shows a high accuracy when compared with manual labelling results. Accordingly, the mean regional asymmetries in the UK Biobank might partly reflect left–right differences present in the reference dataset used to construct the atlas. However, our study was focused primarily on comparing relative asymmetry between genotypes, at the regional level. The use of an asymmetrical atlas based on healthy individuals had the advantage that regional identification was likely to be accurate for structures that are asymmetrical in the general population, while taking hemisphere-specific information into account.

35We also made use of continuous variables as covariates in heritability estimation and genome-wide association analysis (below), which were: age when attended assessment centre (UK Biobank fields 21003-2.0), nonlinear age, that is (age − mean_age)2, the first ten genetic principal components capturing population genetic diversity (fields 22009-0.1 to 22009-0.10) (or the first 40 principal components in the replication dataset with higher diversity of ancestry), scanner position parameters (X, Y and Z position: fields 25756-2.0, 25757-2.0 and 25758-2.0), T1 signal-to-noise ratio (field 25734-2.0) and T1 contrast-to-noise ratio (field 25735-2.0), plus categorical covariates which were: assessment centre (field 54-2.0), genotype measurement batch (field 22000-0.0) and sex (field 31-0.0).

36From the primary dataset, 9,516,074 autosomal variants with minor allele frequencies >1%, INFO score >0.7 and Hardy–Weinberg equilibrium P > 1 × 10−7 were used to build a genetic relationship matrix using GCTA (version 1.93.0beta). Specifically for analyses using GCTA, we further excluded one random participant from each pair having a kinship coefficient higher than 0.025 based on the calculated genetic relationship matrix (as this analysis is especially sensitive to higher levels of relatedness), resulting in 30,315 participants. Genome-based restricted maximum likelihood (GREML) analyses were performed to estimate the SNP-based heritability for each AI, controlling for the above-mentioned covariates, and applying FDR 0.05 across the 73 AIs to define significantly heritable AIs. Bivariate GREML analysis was used to estimate genetic correlations between pairs of AIs, separately for cortical surface area, cortical thickness and subcortical volume AIs, with FDR correction at 0.05 for multiple testing.

37In mvGWAS, a single association test is performed for each SNP in relation to multiple traits simultaneously. We used MetaPhat software to perform mvGWAS analysis across asymmetries for cortical surface area, cortical thickness and subcortical volume, including only the 42 AIs that had shown significant SNP-based heritability. MetaPhat performs meta-canonical correlation analysis, and uses univariate GWAS summary statistics as input from each separate AI, which were derived under an additive genetic model while controlling for the above-mentioned covariates, using BGENIE software (v1.2). Thus our mvGWAS tested effectively for association with 42 traits. This approach estimates the linear combination of traits that is maximally associated with genotype, which can differ for each SNP, while maintaining a correct false-positive rate. A total of 9,803,522 SNPs (see further above) were used for mvGWAS, spanning all autosomes and chromosome X. Statistically significant SNPs were considered as those with P < 5 × 10−8 in mvGWAS, which is a widely used threshold to account for multiple testing over the whole genome, in the context of LD in European-descent populations,.

38MetaPhat also uses systematic criteria to define central traits which make the greatest contributions to significant multivariate associations, on the basis of an iterative process to optimize multivariate model properties with reference to canonical correlation analysis P values and the Bayesian information criterion. For the lead SNPs at genome-wide significant loci (see below for how these were defined), we also performed post hoc analysis in which we examined their separate left and right hemispheric associations, using traits corresponding to the central AIs that were involved in the multivariate associations (Supplementary Table 11), again using BGENIE, an additive genetic model and covariates as described above.

39As a sensitivity analysis, we re-ran the mvGWAS after excluding from the primary dataset 886 participants who had lifetime diagnoses of neurological conditions that could potentially disrupt brain structure (Supplementary Table 20). The significant mvGWAS loci were minimally affected by this exclusion (Supplementary Fig. 9).

40FUMA (version v1.3.6), an online platform for functional annotation of GWAS results, was applied to the results from mvGWAS. A multi-step process, using default parameters, was used to identify distinct, significantly associated genomic loci, and independent lead SNPs within those loci. Briefly, on the basis of the pre-calculated LD structure from the 1000 Genomes European reference panel, SNPs with genome-wide significant mvGWAS P values <5 × 10−8 that had LD r2 < 0.6 with any others were identified. For each of these SNPs, other SNPs that had r2 ≥ 0.6 with them were included for further annotation (see below), and independent lead SNPs were also defined among them as having low LD (r2 < 0.1) with any others. If LD blocks of significant SNPs are located within 250 kb of each other (default parameter), they are merged into one genomic locus. Therefore, some genomic loci could include one or more independent lead SNPs (Table 1). The major histocompatibility complex region on chromosome 6 was excluded from this process by default, because of its especially complex and long-range LD structure.

41Functional annotations were applied by matching chromosome location, base-pair position, reference and alternate alleles to databases containing known functional annotations, which were ANNOVAR categories, Combined Annotation-Dependent Depletion scores, RegulomeDB scores and chromatin state,:ANNOVAR categories identify SNPs on the basis of their locations with respect to genes, such as exonic, intronic and intergenic, using Ensembl gene definitions.Combined Annotation-Dependent Depletion scores predict deleteriousness, with scores higher than 12.37 suggesting potential pathogenicity.RegulomeDB scores integrate regulatory information from eQTL and chromatin marks, and range from 1a to 7, with lower scores representing more importance for regulatory function.Chromatin states show the accessibility of genomic regions, and were labelled by 15 categorical states on the basis of five chromatin marks for 127 epigenomes in the Roadmap Epigenomics Project, which were H3K4me3, H3K4me1, H3K36me3, H3K27me3 and H3K9me3. For each SNP, FUMA calculated the minimum chromatin state across 127 tissue/cell-type in the Roadmap Epigenomics Project. Categories 1–7 are considered as open chromatin states.

42ANNOVAR categories identify SNPs on the basis of their locations with respect to genes, such as exonic, intronic and intergenic, using Ensembl gene definitions.

43Combined Annotation-Dependent Depletion scores predict deleteriousness, with scores higher than 12.37 suggesting potential pathogenicity.

44RegulomeDB scores integrate regulatory information from eQTL and chromatin marks, and range from 1a to 7, with lower scores representing more importance for regulatory function.

45Chromatin states show the accessibility of genomic regions, and were labelled by 15 categorical states on the basis of five chromatin marks for 127 epigenomes in the Roadmap Epigenomics Project, which were H3K4me3, H3K4me1, H3K36me3, H3K27me3 and H3K9me3. For each SNP, FUMA calculated the minimum chromatin state across 127 tissue/cell-type in the Roadmap Epigenomics Project. Categories 1–7 are considered as open chromatin states.

46We also used FUMA to annotate independent significant SNPs and their candidate SNPs according to previously reported phenotype associations (P < 5 × 10−5) in the National Human Genome Research Institute–European Bioinformatics Institute catalogue.

47For a significant mvGWAS association in the major histocompatibility complex region (Table 1), we took the most significant individual SNP, rs7781 (P = 1.62 × 10−10), as the single lead SNP to represent this locus, and annotated it manually.

48SNP-to-gene mapping at significant mvGWAS loci was performed using the default FUMA processes for these three strategies:Positional mapping was used to map SNPs to protein-coding genes on the basis of physical distance (within 10 kb) in the human reference assembly (GRCh37/hg19).eQTL mapping was used to annotate SNPs to genes (that is, when SNP genotypes are associated with variation in gene mRNA expression levels). eQTL mapping was carried out in relation to genes up to 1 Mb away on the basis of four brain-expression data repositories: PsychENCORE, CommonMind Consortium, BRAINEAC and GTEx v8 Brain. FUMA applied a FDR of 0.05 within each analysis to identify significant eQTL associations.Chromatin interaction mapping was performed to map SNPs to genes on the basis of seven brain-related Hi-C chromatin conformation capture datasets: PsychENCORE EP link (one way), PsychENCORE promoter anchored loops, HiC adult cortex, HiC foetal cortex, HiC (GSE87112) dorsolateral prefrontal cortex, HiC (GSE87112) hippocampus and HiC (GSE87112) neural progenitor cell. We further selected only those genes for which one or both regions involved in the chromatin interaction overlapped with a predicted enhancer or promoter region (250 bp up- and 500 bp downstream of the transcription start site) in any of the brain-related repositories from the Roadmap Epigenomics Project, that is, E053 (neurospheres) cortex, E054 (neurospheres) ganglion eminence, E067 (brain) angular gyrus, E068 (brain) anterior caudate, E069 (brain) cingulate gyrus, E070 (brain) germinal matrix, E071 (brain) hippocampus middle, E072 (brain) inferior temporal lobe, E073 (brain) dorsolateral prefrontal cortex, E074 (brain) substantia nigra, E081 (brain) foetal brain male, E082 (brain) foetal brain female, E003 embryonic stem (ES) H1 cells, E008 ES H9 cells, E007 (ES-derived) H1 derived neuronal progenitor cultured cells, E009 (ES-derived) H9 derived neuronal progenitor cultured cells and E010 (ES-derived) H9 derived neuron cultured cells. A FDR of 1 × 10−6 was applied to identify significant interactions (default parameter), separately for each analysis.

49Positional mapping was used to map SNPs to protein-coding genes on the basis of physical distance (within 10 kb) in the human reference assembly (GRCh37/hg19).

50eQTL mapping was used to annotate SNPs to genes (that is, when SNP genotypes are associated with variation in gene mRNA expression levels). eQTL mapping was carried out in relation to genes up to 1 Mb away on the basis of four brain-expression data repositories: PsychENCORE, CommonMind Consortium, BRAINEAC and GTEx v8 Brain. FUMA applied a FDR of 0.05 within each analysis to identify significant eQTL associations.

51Chromatin interaction mapping was performed to map SNPs to genes on the basis of seven brain-related Hi-C chromatin conformation capture datasets: PsychENCORE EP link (one way), PsychENCORE promoter anchored loops, HiC adult cortex, HiC foetal cortex, HiC (GSE87112) dorsolateral prefrontal cortex, HiC (GSE87112) hippocampus and HiC (GSE87112) neural progenitor cell. We further selected only those genes for which one or both regions involved in the chromatin interaction overlapped with a predicted enhancer or promoter region (250 bp up- and 500 bp downstream of the transcription start site) in any of the brain-related repositories from the Roadmap Epigenomics Project, that is, E053 (neurospheres) cortex, E054 (neurospheres) ganglion eminence, E067 (brain) angular gyrus, E068 (brain) anterior caudate, E069 (brain) cingulate gyrus, E070 (brain) germinal matrix, E071 (brain) hippocampus middle, E072 (brain) inferior temporal lobe, E073 (brain) dorsolateral prefrontal cortex, E074 (brain) substantia nigra, E081 (brain) foetal brain male, E082 (brain) foetal brain female, E003 embryonic stem (ES) H1 cells, E008 ES H9 cells, E007 (ES-derived) H1 derived neuronal progenitor cultured cells, E009 (ES-derived) H9 derived neuronal progenitor cultured cells and E010 (ES-derived) H9 derived neuron cultured cells. A FDR of 1 × 10−6 was applied to identify significant interactions (default parameter), separately for each analysis.

52Genome-wide gene-based association analysis was performed using mvGWAS summary statistics as input into MAGMA (v1.08), using default parameters implemented in FUMA (SNP-wide mean model). This process examines the joint association signals of all SNPs within a given gene (including 50 kb upstream to 50 kb downstream of the gene), while considering the LD between the SNPs. SNPs were mapped to 20,146 protein-coding genes on the basis of National Center for Biotechnology Information build 37.3 gene definitions, and each gene was represented by at least one SNP. We applied a Bonferroni correction for the number of tested genes (P < 0.05/20,146).

53We used MAGMA, again with default settings as implemented in FUMA, to test for enrichment of association within predefined gene sets. This process tests whether gene-based P values among all 20,146 genes are lower for those genes within pre-defined functional sets than the rest of the genes in the genome, while correcting for other gene properties such as the number of SNPs. A total of 7,343 gene sets, defined according to Gene Ontology biological processes, were tested from MSigDB version 7.0. In the main text we report the gene sets with P values that met Bonferroni correction for multiple testing (P < 0.05/7,343).

54In addition, we used the list of single closest genes to the 27 lead SNPs arising from mvGWAS (Table 1) as input for gene set enrichment analysis, using the same 7,343 Gene Ontology biological process gene sets, but now on the basis of the hypergeometric test as implemented in GENE2FUNC of FUMA, which is appropriate for gene lists.

55Finally, we used the CELL TYPE function (as implemented within FUMA) to test whether lower gene-based association P values for brain asymmetry were associated with differential expression levels across cell types, using Bonferroni correction within each separate analysis with respect to each cell-type expression dataset included in FUMA.

56We used the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://string-db.org) for protein network analysis, using as input the names of 57 genes identified through gene-based association analysis, as described above. The STRING dataset includes protein–protein interaction information from numerous sources, including experimental data, publications and computational prediction methods. Only links with medium confidence or higher (confidence score >0.4; default parameter) were retained.

57Using the gene-based association P values for all 20,146 genes genome wide, we used MAGMA (default settings as implemented in FUMA) to examine whether lower gene-based P values tended to be found for genes showing relatively higher expression in BrainSpan gene expression data from any particular ages compared with all other ages, separately for 29 different age groups ranging from 8 postconceptional weeks to 40 years old, and 11 defined developmental stages from early prenatal to middle adulthood. We corrected for multiple testing through a FDR of 0.05 (separately for the two analyses).

58We applied the iSECA toolbox that can test for genetic overlap on the basis of per-SNP association P values only (mvGWAS does not produce univariate beta coefficient effect size estimates that can be used in standard genetic correlation analysis). We tested for genetic overlap in relation to traits previously reported to associate phenotypically with different aspects of brain structural asymmetry (Main), using GWAS P values from previously published, large-scale studies: educational attainment (n = 1,131,881), handedness (n = 331,037), intelligence (n = 269,867), autism spectrum disorder (n = 46,350), attention deficit hyperactivity disorder (n = 55,374) and schizophrenia (n = 82,315). We also tested for genetic overlap in relation to brain intracranial volume (n = 32,438). After LD-based filtering and clumping using default parameters, iSECA tests for pleiotropy between two sets of GWAS results using an exact binomial statistical test at each of 12 P value levels: P ≤ (0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1). The analysis compares the expected and observed overlap in the subsets of SNPs at these levels from two GWAS (144 combinations in total). In other words, iSECA iterates through each of the 12 P value levels and counts the number of overlapping variants between two GWAS at each P value threshold, and compares that number with the number expected under the null hypothesis of no genetic overlap, using the exact binomial test. iSECA then counts up the number of comparisons with evidence of overlap at a nominally significant level of P ≤ 0.05. To assess the significance level of overlap, we generated 1000 datasets through permutations (default parameter), which contained all the possible combinations for a pair of traits, and determined whether the number of levels with nominally significant genetic overlap was significantly more than expected by chance. Finally, Bonferroni correction <0.05 was applied for multiple testing of seven traits. Additionally, iSECA generated Q–Q plots for asymmetry mvGWAS P values conditioned on the other trait P values (for example, P ≤ 0.1, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0) to visualize whether there is an excess of pleiotropic SNPs, which should be visible as a leftward shift of the curve as the P value threshold becomes tighter (Supplementary Fig. 8).

59Further information on research design is available in the Nature Research Reporting Summary linked to this article.