The genetic architecture of human cerebellar morphology supports a key role for the cerebellum in human evolution and psychopathology

Data-driven decomposition of cerebellar grey matter maps reveals highly reproducible structural covariance patterns

Since traditional atlases of the cerebellar cortex based on gross anatomical landmarks (i.e., lobules) only partially overlap with more recent functional parcellations of the cerebellum^17–25, we first used a data-driven approach (non-negative matrix factorization, NNMF^26–28) to parcellate voxel-based morphometry (VBM) based maps of cerebellar grey matter volume (containing 147,121 1*1*1 mm voxels) from 28,212 participants into a smaller number of structural covariance patterns (SCPs), i.e., cerebellar sub-regions where voxel-volumes co-vary consistently across individuals (see Online Methods for details about the study sample, MRI processing and quality control procedures). Like other dimensionality reduction methods, such as principal component analysis (PCA) or independent component analysis (ICA), NNMF decomposes the input matrix (here 147k voxels * 28k participants) into two lower rank matrices (voxel weights: W, and participant weights: H), so that the product of W and H approximates the original input data. The defining feature of NNMF is that it requires these lower rank matrices to be non-negative, which has previously been shown to result in sparse, reproducible and easily interpretable parcellations of high-dimensional brain imaging data^26–30. In essence, when applied to voxel-wise indices of grey matter volume, NNMF yields distinct maps of voxels that show similar patterns of volume-variation across individuals, commonly referred to as “structural covariance patterns”^26,27,29. For each structural covariance pattern (common across all participants), NNMF also provides individual subject weighs expressing the degree to which these patterns are expressed, i.e., reflecting individual variation in regional cerebellar volumes.

In the current study NNMF yielded highly reproducible cerebellar structural covariance patterns (across split-half datasets) for model orders (i.e., number of components/patterns specified) ranging from 2 to 100 (see Supplementary Figures 2 for summary maps of all tested model orders). After observing that the improved fit to the original data seen with higher model orders tended to level off between 15 and 30 components (indicating that the intrinsic dimensionality of the data might have been reached (see Supplementary Figure 3), we decided on a model order of 23 based on its good split-half reproducibility (see Fig. 1B, and Online Online Methods for details regarding model order selection). Importantly, when subject weights were summed across the 23 structural covariance maps, the resulting values corelated tightly (r = 0.9995) with estimates of total cerebellar grey matter volume (see also Supplementary Figure 4), demonstrating that our data-driven decomposition preserved inter-individual variation in cerebellar volume.

• Download figure
• Open in new tab

Figure 1: Data driven decomposition of cerebellar grey matter maps yields highly reproducible and moderately heritable structural covariance patterns (SCPs).

A: Binarized winner-takes-all map for the 23-component solution based on data-driven decomposition of cerebellar grey matter maps from 28,212 participants. Note that empirically derived boundaries between cerebellar regions only partially follow traditional lobular borders (marked with dotted black lines); B: Five distinct components overlapping cerebellar Crus I derived from the split-half reliability analyses. While one of these SCPs emerged as bilateral in Split 2, the remaining four SCPs were almost identical, despite being derived from two independent samples; C: Narrow sense (SNP-based) heritability of the 23 SCPs (see Supplementary Data 5 for numerical values). D: Hierarchical clustering of the 23 SCPs (derived from the full sample decomposition) based on their pairwise genetic correlations revealed a primary division between the anterior and posterior cerebellum, with additional separations between medial and lateral regions. The full genetic correlation matrix can be found in Supplementary Data 6.

Of note, our data-driven decomposition differed markedly from the standard cerebellar atlas based on gross anatomical features, shown as dotted lines in Fig 1A (see Supplementary Figure 3 for all 23 components and Supplementary Data 1 for quantification of overlap between NNMF-derived structural covariance patterns (SCPs) and standard cerebellar anatomical regions, i.e. lobules). For instance, five distinct SCPs overlapped Crus I of cerebellar lobule VII (shown in Fig 1B), an anatomical region which already started to split into separate components at a model order of three. We further observed only partial overlap with task-based functional parcellations of the cerebellar cortex (Supplementary Data 2). While some SCPs clearly overlapped cerebellar regions previously associated with hand movements, eye-movements/saccades or autobiographical recall, other data-driven SCPs overlapped multiple functionally defined cerebellar regions (Supplementary Data 2).

Cerebellar structural covariance patterns are heritable and reveals a distinct anterior-posterior pattern based on their bivariate genetic correlations

After removal of one of each genetically related pair of individuals (n = 910), 27,302 participants remained for the genetic analyses. In addition to the 23 regional cerebellar structural covariance patterns of primary interest, we also included total cerebellar volume, estimated total intracranial volume and 9 cerebral brain phenotypes to serve as covariates and/or comparison phenotypes. Prior to all genetic analyses, morphological features were adjusted for effects of scanner site, sex, age, estimated total intracranial volume, 40 genetic population components, genetic analysis batch and a quantitative structural MRI quality index (the Euler number³¹) using general additive models, before finally being rank-order normalized (see Online Methods for details).

To validate our analysis approach, we computed genetic correlations (using LD-score regression, LDSC³²) between univariate GWAS results on the comparison features (see Manhattan- and QQ-plots in Supplementary Figure 5) and previously published neuroimaging GWAS studies on these brain phenotypes. Results showed a mean r_g of .90 (range: .80-.99, see Supplementary Data 3). For the 23 cerebellar morphological features, we used the univariate GWAS summary statistics (see Supplementary Figures 6-7 for Manhattan- and QQ-plots) to compute genetic correlations between discovery (n = 27,302) and replication (n = 11,264) samples using LDSC³². These genetic correlations were high (mean rG: .92; range: .83-1), indicating reliable genetic signals (see Supplementary Data 4).

Genetic complex trait analysis (GCTA³³) revealed SNP-based heritability estimates (h²) ranging from .33 to .44 (Figure 1C and Supplementary Data 5). Analyses of total cerebellar volume (h² = .35), estimated total intracranial volume (h² = .41) and the 9 cerebral comparison phenotypes (h² range: .26 to .45) gave a similar range of heritability estimates (Supplementary Data 5).

Hierarchical clustering of cerebellar features based on their bivariate genetic correlation matrix (GCTA-based rG ranging from 0.35 to 0.98, Supplementary Data 8) revealed a primary anterior-posterior division running along the horizontal fissure separating Crus I and Crus II, with secondary divisions grouping features into more lateral or medial, as well as more anterior and posterior features within the major anterior and posterior regions (Fig 1D). Of note, the primary division along the horizontal fissure was also evident from the (genetically naïve) two-component NNMF decomposition, while medial-to-lateral divisions already began to emerge with a model order of three (see Supplementary Figure 1).

In order to examine whether this phenotypic and genetic correlation structure was also reflected in regional cerebellar gene expression patterns, we used the abagen toolbox³⁴ to extract Allen Human Brain Atlas³⁵ gene expression profiles for 22 of the 23 morphological features and computed their bivariate Pearson correlations (across 15,631 genes; Supplementary Data 7) and hierarchical clustering solution (Supplementary Figure 8). The anterior-posterior boundary across the horizontal fissure was also evident in the gene expression data, which in addition highlighted distinct gene expression patterns for the posterior midline (grouped together with the horizontal fissure), as well as for the most lateral regions of the cerebellar cortex.

Multivariate GWAS reveals 351 genetic loci associated with cerebellar morphology

Figure 2A shows the main results for the multivariate GWAS across the 23 cerebellar structural covariance patterns applying MOSTest³⁶. We observed 35,098 genome-wide significant (GWS) SNPs, which were mapped by FUMA³⁷ to a total of 51,803 candidate SNPs by adding reference panel SNPs in high LD (r >.6) with GWS SNPs. The 51,803 candidate SNPs (see Supplementary Data 8) were represented by 1,936 independent significant SNPs and 560 lead SNPs in 351 genomic loci (see Supplementary Data 10). The LDSC intercept of the MOSTest summary statistics was 1.0352 (i.e., indicative of no or minimal genetic inflation), while the QQ-plot of results based on permuted data (under the null hypothesis) confirmed the validity of the MOSTest analytical approach (Supplementary Figure 9).

• Download figure
• Open in new tab

Figure 2: Multivariate GWAS analysis of the 23 cerebellar structural covariance patterns revealed 351 independent genome-wide significant (GWS) loci.

A: The upper half of the Miami plot shows the main results from the multivariate analysis. The lower half displays results from a series of 23 univariate analyses (corrected for multiple comparisons using the standard min-P approach), as well as results from a univariate analysis of total cerebellar grey matter (marked in orange). B-D: Euler diagrams showing the relative numbers of - and overlaps between - candidate SNPs (in thousands) mapped by the three analysis approaches employed in the current study (B), the current and four recent studies reporting genetic associations with cerebellar morphology (C), as well as results from multivariate GWASs on hippocampal and cerebrocortical morphology (D). For full results on overlap between all cerebellar candidate SNPs and other brain phenotypes, see Supplementary Data 8.

Annotation of all candidate SNPs using ANNOVAR³⁸ as implemented in FUMA³⁵ revealed that the majority of candidate SNPs were intronic (57.8%) or intergenic (38.3%). While only 0.7% were exonic, about 81% of the candidate SNPs were assigned minimal chromatin states between 1 and 7 (i.e., open chromatin states), implying effects on active transcription³⁷ (see Supplementary Data 8-9).

We evaluated the robustness of these multivariate results using a multivariate replication procedure established in Loughnan et al³⁹, which computes a composite score from the mass-univariate z-statistics (i.e., applying multivariate weights from the discovery sample to the replication sample input data) and then tests for associations between this composite score and genotypes in the replication sample (for mathematical formulation see Loughnan et al³⁹). Results showed that 97% of the 339 loci lead SNPs present in both samples (i.e., 94% of the 351 reported loci lead SNPs from the discovery sample) replicated at a nominal significance threshold of p < .05 (Supplementary Figure 10 and Supplementary Data 10), and that 74% remained significant after Bonferroni correction (Supplementary Data 10). Moreover, 99% of loci lead SNPs showed the same effect direction across discovery and replication samples (Supplementary Figure 9 and Supplementary Data 10).

In addition, we assessed the robustness of the multivariate patterns by computing bivariate correlations between feature z-score vectors assigned to the discovery sample lead SNPs in an independent multivariate GWAS (MOSTest) performed on the replication sample. These correlations (restricted to the 339 loci lead SNPs present in both samples) were relatively high (mean r: .70, see Supplementary Data 10 and Supplementary Figure 11). Figure 3 and Supplementary Figure 11 also give some examples of discovery and replication sample multivariate patterns projected back onto the cerebellar cortex.

• Download figure
• Open in new tab

Figure 3: Loci lead SNPs show spatially heterogeneous and replicable effects across the cerebellar cortex.

The 351 loci lead SNPs identified by MOSTest are plotted as a function of main overall effect across all cerebellar features (x-axis: mean Z-score) and most extreme effect for any single cerebellar feature (y-axis: most extreme Z-score across features), and colour coded by SNP discovery method. The cerebellar flat-maps show discovery (left) and replication (right) sample regional distributions of Z-scores (color-scale range from - 10 to 10) for a few selected lead SNPs (see Supplementary Data 10 for individual feature Z-scores for all 351 discovery sample loci lead SNPs). SNP rs7877685 was only present in the discovery sample.

To compare our main multivariate MOSTest results to univariate approaches, the lower part of the Miami plot in Figure 2A, and Figure 2B, displays results from a set of univariate GWASs on the cerebellar morphological features (which yielded 8370 candidate SNPs and 57 genomic loci, corrected for multiple comparisons using the min-P approach^40,41, see also Supplementary Figures 6-7 for univariate Manhattan and QQ-plots), as well as the 4044 candidate SNPs and 10 significant loci resulting from the univariate GWAS on total cerebellar grey matter volume (marked in orange). 52 of the 57 loci identified in the univariate analyses of regional cerebellar features overlapped 55 of the 351 loci identified using the multivariate method, while 9 out of 10 loci identified in the univariate analysis of total cerebellar volume were also identified in the multivariate analysis (the slightly mismatching numbers are due to loci of unequal lengths, causing some larger loci to overlap with several smaller loci). Thus, our multivariate analysis of regional cerebellar features increased the locus yield ∼35-fold relative to analyzing total cerebellar volume and ∼6-fold relative to performing a set of univariate analyses on the same regional features.

We next compared the current findings with previously reported genetic loci for cerebellar morphology by extracting summary statistics from two recent GWAS studies using the UKBB sample (n = 19k⁴² and 33k⁴³) that included regional cerebellar volumes among the full set of analyzed brain imaging-derived phenotypes (101⁴² and 3,144⁴³, respectively), as well as two recent GWASs on total cerebellar volume (n = 33k¹³ and 27k¹⁴, respectively). Candidate SNPs and independent GWS loci were identified in FUMA using the same settings as for our primary analyses and employing a liberal p-value threshold of 5e-8 (i.e., not correcting for the total number of brain imaging features analyzed). Results are displayed in Figure 2c and Supplementary Data 8 and 10. In brief, we found that 19,527 of the 51,803 candidate SNP (i.e., 36%) and 123 of the 351 identified genomic loci (i.e., 35%) overlapped with candidate SNPs and loci extracted from these three previous studies. Thus, 228 of the 351 (i.e., 65%) genetic loci reported here are novel to the literature on of cerebellar morphology genetics (see Supplementary Data 10.

Overlap of cerebellar candidate SNPs and genetic loci with results from recent multivariate analyses of hippocampal and cerebrocortical morphology are displayed in Figure 2D and Supplementary Data 8 and 10 (final columns). Of note, we found that 32 and 29 percent of the candidate SNPs discovered here for the cerebellum overlapped with candidate SNPs for vertex-wise cerebrocortical surface area and thickness⁴⁴, respectively, while 11.4 percent overlapped with candidate SNPs found for hippopcampal subregions⁴⁵. 95 of the 351 genetic loci overlapped loci linked to the other multivariate brain phenotypes (Supplementary Data 10). Thus, 64% of the candidate SNPs and 73% of genetic loci appeared to be specifically associated with cerebellar morphology.

Significant genetic variants show heterogeneous effects across the cerebellar cortex, influencing both regional and total volumes

A major advantage of our multivariate analysis approach is its sensitivity to both highly localized and more generally distributed effects of SNPs on cerebellar morphology. This is illustrated in Figure 3, which displays the 351 loci lead SNPs as a function of both the most extreme individual Z-score across all cerebellar features (e.g, analogous to the strongest “local” effect) and of the mean Z-score across these features (i.e. analogous to the main effect on overall cerebellar volume).

As can be seen, some loci lead SNPs (e.g. rs13107325; rs76934732) show pronounced positive or negative mean z-scores, indicating a relatively consistent direction of effect across cerebellar features. See also inset figures displaying feature Z-scores projected back onto the cerebellar cortex. Many of these SNPs also emerged in the univariate analysis of total cerebellar volume and were recently reported in GWASs on total cerebellar volume^13,14.

Many other loci lead SNPs, however, show strong “local” signals with opposite effect directions across features, yielding very weak global signals (e.g. rs117332043; rs78777685). Thus, while several of the most significant SNPs in this category have previously been reported in GWASs including local cerebellar morphological features^42,43, they did not emerge from analyses of total cerebellar volume, neither in the current nor in previous studies.

Finally, our multivariate MOSTest approach is also sensitive to SNPs displaying weaker effects distributed across several cerebellar features, often with opposite effect directions (e.g. rs12464825; rs2388334). For this category of SNPs, neither of the univariate methods have sufficient power at the current sample size. In contrast, these two example SNPs robustly emerged from the multivariate analysis (discovery sample p-values < 1e-56; replication sample p-values <1e-15).

Genetic variants associated with cerebellar morphology are enriched for evolutionary recent mutations in the human genome

We next mapped the evolutionary age of cerebellar lead SNPs (thresholded for linkage disequilibrium at r2 <0.1) by merging them with a recently published dataset on dated mutations in the human genome⁴⁶. Following the analysis procedure established by Libedinsky et al.⁴⁷, we plotted the histogram of dated lead SNPs over the last 2 million years in bins of 20.000 years (Figure 4A and Supplementary Data 11), and tested for positive or negative enrichment by comparing them to empirical null distributions derived from 10,000 randomly drawn and equally sized sets of all SNPs in the full human genome dating dataset (after matching them to cerebellum-associated lead SNPs in terms of minor allele frequency; see Online Methods for details).

• Download figure
• Open in new tab

Figure 4: Lead SNPs associated with cerebellar morphology are enriched for evolutionary recent mutations in the human genome.

A: Histogram of estimated SNP age (ranging from 0 to 2 million years, in bins on 20,000 years) for 548 independent lead SNPS associated with cerebellar morphology. The solid black and dotted lines denote the mean and upper 95th confidence interval (Bonferroni corrected across 100 time-bins) derived from a null model constructed from 10,000 of equally sized sets of SNPs randomly drawn from the Human Genome Dating Atlas of Variant Age (after matching these to cerebellar lead SNPs in terms of minor allele frequencies). Red bars denote time-bins of significant positive enrichment. B: Histogram showing comparative evolutionary enrichment effects for lead SNPs identified in a multivariate GWASs of cerebrocortical area (for cerebrocortical thickness, see Supplementary Figure 12), See Supplementary Data 11-12 for full numerical results.

Results revealed positive enrichment for cerebellar lead SNPs in the time bin ranging from 20-40,000 years ago, i.e., overlapping the Upper Paleolithic, a period characterized by rapid cultural evolution and the first evidence of several uniquely human behaviors (often referred to as behavioral modernity), such as the recording of information onto objects⁴⁸.

For comparison, Figure 4B also shows results from an analysis of lead SNPs identified in a previous multivariate GWAS study of regional cerebrocortical surface area ⁴⁴, while Supplementary Figure 12 shows the corresponding histogram for regional cerebrocortical thickness. As can be seen, both these cerebrocortical phenotypes also showed significant enrichment in overlapping time bins (i.e., 20-60,000 years ago). For full numerical results in Supplementary Data 12. Given the relatively low number of independent lead SNPs (range: 548-862 SNPs across phenotypes), we also ran validation analyses using all independent significant SNPs (LD-thesholded at r2 < 0.6, range: 1574-2883 SNPs), which yielded very similar results (see Supplementary Figure 13 and Supplementary Table 13).

Genes associated with cerebellar morphology show selective expression in cerebellar and prenatal brain tissue, as well as enrichment for genes linked to human accelerated regions

To functionally characterize the multivariate GWAS signal, we mapped the full set of GWAS p-values to 19,329 protein coding genes using MAGMA⁴⁹ and used the resulting gene-level p-values to test for 1) GWS genes, 2) gene expression in brain tissue; and 3) enrichment for genes linked to human accelerated regions (HARs), i.e. sections of DNA that have remained relatively conserved throughout mammalian evolution, before being subject to a burst of changes in humans since divergence of humans from chimpanzees^50,51. These analyses yielded a total of 534 GWS genes (i.e., 2.78% of all protein coding genes, see Figure 5A and Supplementary Data 14). Using the full set of 19,329 gene-level p-values in MAGMA gene property analyses revealed significant and specific gene expression in cerebellar and prenatal brain tissue (Figure 5B-C and Supplementary Data 15-16), with the selective cerebellar expression seen in two independent datasets (Allen Human Brain Atlas³⁵ and The Genotype-Tissue Expression (GTEx) Project).

• Download figure
• Open in new tab

Figure 5: Gene mapping reveals selective enrichment for cerebellar and prenatal brain tissue, as well as for genes linked to human accelerated regions.

A: MAGMA mapped the full range of SNPs from the multivariate GWAS to 19.329 protein-coding genes, 534 of which were genome-wide significant (Bonferroni corrected threshold marked with red line), with the 30 most significant genes labeled in the gene-based Manhattan plot. B-C: MAGMA gene-property analyses revealed selective brain expression in the cerebellum and during prenatal developmental stages. D: MAGMA gene-set analysis revealed significant enrichment for sets of genes previously linked to human accelerated regions (HARs). Figure 4d also shows results from comparative analyses using summary statistics from other multivariate GWASs of MRI-based brain, as well as significant results from statistical tests comparing beta-weights for HAR-linked genes across multivariate brain features. Horizontal red lines mark the Bonferroni-corrected significance threshold for each subplot. See Supplementary Data 14-17 for full results.

The MAGMA gene set analysis of HAR-linked genes revealed significant enrichment (p = 7.09e-08) for genes associated with cerebellar morphology (Figure 5D, Supplementary Data 17). Of note, running this same HAR gene set analysis on summary statistics from recent multivariate GWAS studies on cerebrocortical⁴⁴ or hippocampal⁴⁵ morphological features yielded similar (for cortical features) or significantly weaker (hippocampal features) enrichment effects (Figure 5d, see Online Methods for processing pipeline).

Genes linked to human cerebellar morphology show enrichment for gene sets linked to selective cerebellar gene expression, altered cerebellar morphology in mouse models, human clinical/anthropometric traits, as well as specific biological processes and pathways

In line with the continuous brain tissue gene expression results described above, we also observed significant and relatively selective enrichment for smaller curated sets of genes previously found to be highly and selectively expressed in mouse (Figure 6A, Supplementary Data 18) and human (Figure 6C, Supplementary Data 20) cerebellar brain tissue, as well as for sets of genes previously shown to affect cerebellar morphology in mouse gene perturbation experiments (Figure 6B, Supplementary Data 19) and various clinical conditions and anthropometric traits in humans (Figure 6D, Supplementary Data 21). The most significant gene ontology and curated gene sets from the MSigDB^52,53 database were related to brain development (e.g., neurogenesis, axon guidance and neuron differentiation, Figure 6E, Supplementary Data 22), and highlighted the reelin signaling pathway (Figure 6F, Supplementary Data 23).

• Download figure
• Open in new tab

Figure 6: Gene mapping reveals selective enrichment across brain tissues and curated gene sets.

MAGMA gene-set analyses revealed significant enrichment for sets of genes previously linked to preferential expression in human (A) and mouse (C) cerebellar tissue, as well as effects on cerebellar morphology (and other brain phenotypes) in mouse gene perturbation experiments (B) and human clinical disorders and phenotypes (D). Significant gene ontology terms were related to neural development (E), while curated gene sets highlighted the Reelin signaling pathway (F). Across all subplots the x-axis shows the -log10 p-value, the y-axis marks the number of gene sets, and the top 10 most significant gene-sets are labelled. Red lines mark the Bonferroni-corrected significance threshold for each subplot. See Supplementary Data 18-23 for full results.

As can be seen in Figure 5A (and Supplementary Data 14), RELN (encoding the protein Reelin) was also the most significant gene mapped by MAGMA. See also Supplementary Figure 14 for a regional locus plot showing the 12 lead SNPs mapped to RELN and their associated z-score maps.

Gene mapping reveals sets of plausible causal genes

In addition to the gene-based mapping strategy using all SNPs described above (MAGMA), we also mapped candidate SNPs to genes using two complementary gene mapping strategies: 1) positional mapping of deleterious SNPs (defined as having a CADD-score⁵⁴ > 12.37); and 2) mapping of SNPs previously shown to affect gene expression in cerebellar tissues (i.e., eQTL mapping). Across all three strategies, we mapped a total of 674 unique genes; 531 using MAGMA, 310 using positional and 197 using eQTL mapping (see Figure 7A and Supplementary Data 24). 298 genes were identified by at least two strategies, while 65 genes were mapped by all three strategies. Out of these 674 genes, 61 have previously been associated with cerebellar pathology in humans and/or altered cerebellar morphology in mouse gene perturbation experiments, while 121 have been linked to human accelerated regions (Supplementary Data 24). As can be seen in Figure 7B and Supplementary Data 24, the 674 genes mapped to cerebellar morphology showed some overlap with genes mapped to hippocampal⁴⁵ and cerebrocortical⁴⁴ morphology using the same mapping strategies, but 264 genes (39%) appeared relatively specific to the cerebellum.

• Download figure
• Open in new tab

Figure 7: Mapping and functional characterization of plausible causal genes.

The 674 plausible causal genes mapped using three complementary strategies (A) show partial overlap with genes mapped to hippocampal⁴⁵ and cerebrocortical⁴⁴ morphology (B); significant enrichment for sets of genes linked to human accelerated regions (HAR) (C); selective expression in the cerebellum in mice (D) and humans (E), as well as effects on cerebellar morphology in mouse gene perturbation studies (F) and human disorders and anthropometric phenotypes (G). Across all subplots the x-axis shows the -log10 p-value, the y-axis marks the number of gene sets, and the top 10 most significant gene-sets are labelled. Red lines mark the Bonferroni-corrected significance threshold for each subplot. See Supplementary Data 24-29 for full results.

Results from gene set analyses on the 674 mapped genes (Figure 7, Supplementary Data 25-29) largely mirrored results from the MAGMA analyses described above, but in addition revealed that this set of 674 mapped genes was also enriched for gene sets associated with several complex clinical phenotypes and anthropometric traits in humans, including cognitive ability, neuroticism and schizophrenia (Figure 7F, Supplementary Data 29).

Of note, while the full set of mapped genes showed significant enrichment for sets of genes known to alter cerebellar morphology in mouse mutation or knock-down experiments (Supplementary Data 27), we also note that 613 of the 674 mapped genes have not to our knowledge previously been linked to cerebellar development, anatomy or pathology in mice or humans (Supplementary Data 24), and thus constitute potential targets for future gene perturbation experiments in animal models.

Restricting the above analyses to the 298 genes mapped across at least two strategies did not markedly affect the results (Supplementary Data 25-29, final columns).

Cerebellar morphology shows significant genetic overlap with psychiatric disorders

We finally tested for overlap between the multivariate genetic profile for cerebellar morphology and genetic profiles for five major developmental/psychiatric disorders (attention deficit hyperactivity disorder: ADHD; autism spectrum disorder: ASD; bipolar disorder: BIP; major depressive disorder: MDD; and schizophrenia: SCZ) using conditional/conjunctional FDR analysis⁵⁵. As can be seen in the conditional QQ-plots in Figure 8, these analyses revealed clear patterns of enriched association with the clinical phenotypes when selecting subsets of SNPs with increasingly stronger association with cerebellar morphology (Fig. 8A and figure insets in (C-F).

• Download figure
• Open in new tab

Figure 8: Genetic variants influencing cerebellar morphology overlap with variants associated with five major mental disorders.

Conditional QQ plots (A and figure insets in C-F) show an incremental incidence of association with five mental disorders (leftward deflection) as a function of the significance of association with cerebellar morphology. Manhattan plots show SNPs with significant association with both traits, thresholded at a conjunctional FDR threshold of p > .01 (red dotted line). SCZ: Schizophrenia; BIP: Bipolar Disorder; MDD: Major Depressive Disorder; ADHD: Attention Deficit Hyperactivity Disorder; ASD: Autism Spectrum Disorders.

See also Supplementary Figure 15 for QQ-plots depicting the reverse association, i.e., enriched association with cerebellar morphology when conditioning on the association with developmental/psychiatric disorders. Specific genetic variants jointly influencing the two phenotypes were identified using conjunctional FDR analyses at a conservative statistical threshold of p<.01. Results revealed shared genetic loci with all disorders; namely 48 with SCZ, 22 with BIP, 2 with MDD, 5 with ADHD and 5 with ASD (Figure 6; Supplementary Data 30-34). We mapped lead and candidate SNPs for each of these loci to genes using positional and eQTL mapping and checked for gene overlap across disorders (Supplementary Data 35). Of note, the LRP8 gene (a HAR-linked gene^50,51 encoding a Reelin receptor) emerged from the conjunctional FDR analyses of both BIP and SCZ, thus again highlighting the Reelin signaling pathway.

To complement these multivariate analyses, we also conducted a set of univariate analyses, using LD-score regression³² to calculate the genetic correlations between each individual cerebellar feature and the five developmental/psychiatric disorders.

As can be seen in the first row of Figure 9 (and Supplementary Data 36), genetic correlations with cerebellar morphological features were predominantly negative across diagnoses, indicating that genetic variants associated with a clinical diagnosis tended to also be associated with reduced cerebellar volumes (with 42 out of 115 tested associations showing nominally significant negative correlations; Fig 9, second row). However, all univariate genetic correlations were relatively weak, and only a few negative genetic correlations with BIP and SCZ survived Bonferroni correction for the 23 features tested. When also correcting for the five clinical conditions, only the negative correlation between BIP and cerebellar feature 23 (primarily overlapping vermal lobules VIIIa and VIIIb) remained significant.

• Download figure
• Open in new tab

Figure 9: Univariate genetic correlations between cerebellar morphological features and five major mental disorders are negative.

The top row display unthresholded genetic correlations, while these are filtered at increasingly strict statistical thresholds in the following rows, i.e. nominal p-value < .05 (second row); Bonferroni correction for the 23 cerebellar features tested (third row); and Bonferroni correction for both 23 features and 5 clinical conditions (bottom row). Black dotted lines denote lobular boundaries.

Supplementary Data 37 show genetic correlations between the five disorders and the ten comparison brain phenotypes, as well as total cerebellar volume. In general agreement with the observed pattern for regional cerebellar features, total cerebellar grey matter volume showed nominally significant negative genetic correlations with BIP (rg: -.10; p = 0.0052) and SCZ (rg: -.09; p = 0.0133). Pallidal volume also showed nominally significant negative genetic correlations with these two disorders (BIP: rg: -.10; p = 0.0087; SCZ: rg: -.08; p = 0.0172, while ADHD displayed negative genetic correlations with estimated total intracranial volume (rg: -.15; p = 0.0003) and total cortical surface area (rg:-.14 ; p = 0.0065), as well as a positive genetic correlation with hippocampal volume (rg: .13; p = 0.0164). Finally, MDD showed nominally significant negative genetic correlations with volumes of the hippocampus (rg -.08: p = 0.025) and thalamus (rg: -.09; p = 0.0258). Only the negative genetic correlation between ADHD and estimated total intracranial volume survived Bonferroni correction across the 55 tests performed.

Participants

For our main analyses, T1-weighted MR images, demographic and genetic data from 39,178 UK Biobank participants were accessed using access number 27412. After removing 1043 participants who either had missing genetic data or had withdrawn consent (as of 19.11.2019), data from 38,135 participants remained for the main analysis (age range: 44.6-82.1; mean age: 64.1, 51.9% female). Following quality control procedures (QC, see below), 28,212 UK Biobank participants of European descent remained for the main analyses (age range: 45.1-82.1; mean age: 64.1, 55.1% female), 910 of which were identified as close relatives and removed prior to genetic analyses (see below), leaving a final sample for the primary analyses of 27,302 (age range: 45.1-82.1; mean age: 64.1, 54.9% female). For the replication sample, we accessed a newer release of UK Biobank participants (n = 48,045) and removed the 27,302 participants included in the primary analyses. After running through identical QC procedures as for the main sample (although applied only to the cerebellar features of primary interest and covariates included when analyzing these), the replication sample consisted of 11,264 unrelated UK Biobank participants of European descent (age range: 46.1-83.7; mean age: 66.8, 46.9% female). The UK Biobank was approved by the National Health Service National Research Ethics Service (ref. 11/NW/0382),

Initial MR image processing

MRI data was first processed using the recon-all pipeline in Freesurfer 5.3⁷³, yielding a large number of brain features. Of these, we retained measures of estimated total intracranial volume (eTIV), total cerebro-cortical surface area, average cerebro-cortical thickness as well as the volumes of seven subcortical structures (hippocampus, amygdala, thalamus, pallidum, putamen, caudate nucleus and nucleus accumbens) as cerebral (or global, in the case of eTICV) comparison regions for our main cerebellar analyses. These anatomical features were averaged across hemispheres, yielding a total of ten comparison phenotypes.

Next, the bias-field corrected T1-images from the FreeSurfer analyses were analyzed using the cerebellum-optimized SUIT-toolbox⁷⁴. In brief, SUIT isolates the cerebellum and brain stem from T1-images, segments cropped images into grey and white matter, and normalizes these tissue probability maps to a cerebellum-specific anatomical template. After multiplying the grey matter maps with the Jacobian of the transformation matrix (i.e. preserving information about absolute volumes), we extracted grey matter intensity values overlapping the 28 cerebellar lobular labels in the probabilistic SUIT-atlas. Since cerebellar volumetric indices showed high correlations across the two hemispheres (range: .87-.96; mean: .94), we created mean volumetric measures by averaging across hemispheres. Finally, we also combined the two smallest regions in the SUIT-atlas (Vermis Crus I and II located in the midline, average volumes: 2.9 and 293.4 mm², respectively) to create a new Vermis Crus region (average volume: 296.2 mm²). This procedure reduced the 28 cerebellar lobular volumes to 16 morphological indices. Total cerebellar grey matter volume was defined as the sum of all 28 lobular labels in the SUIT-atlas.

Quality control procedures

After excluding 4,818 UKBB participants of non-European origin, anatomical indices from the remaining 33,317 participants went through an iterative quality control (QC) procedure. First, we excluded 639 subjects with a mean Euler number below -217, indicating poor MRI quality³¹, as well as 12 subjects with missing and 90 subjects with zero values for any of the key cerebellar or cerebral brain measures of interest. Next, we used general additive models (GAM, implemented in the R-package “mgcv”) in order to model the effects of age (estimated as smooth functions for males and females separately, using cubic splines with 10 knots), sex, and scanner site on estimated total intracranial volume (eTIV) and mean thickness of the cerebral cortex. Specifically, we used the following formula:

Adjusted eTIV and mean cortical thickness indices were then created by reconstructing the data using the intercept and residuals from this model (i.e., removing effects of age, sex and scanner), before identifying and rejecting potential outliers, defined as +/− 3 median absolute deviations (MAD)⁷⁵ from the median of these adjusted values. Data from 350 subjects were rejected based on these criteria.

For the remaining cerebral anatomical measures, as well as total cerebellar volume, this procedure was then repeated, with scaled eTIV as an additional predictor. Since previous studies have demonstrated that the relationship between regional brain volumes and intracranial volume is not strictly allometric^76,77, we estimated the effect of eTIV using cubic splines with 10 knots. Specifically, we used the following formula:

In order to be maximally sensitive to outliers in relative cerebellar volumes, we replaced eTIV with total cerebellar volume in the GAM models of cerebellar regions of interest (i.e., SUIT atlas regions), using the following formula:

Adjusted cerebral and cerebellar indices were then created by reconstructing the data using the intercept and residuals from these models (i.e., removing estimated effects of age, sex and eTIV or total cerebellar volume), before rejecting 1792 participants with potential outlier cerebral indices and 2222 participants with potential outlier cerebellar indices (MAD > +/− 3).

This iterative QC procedure resulted in the rejection of 5,105 (i.e., 15.3%) of the original 33,317 datasets, leaving 28,212 datasets for further analysis.

Non-negative matrix factorization

Cerebellar grey matter maps from 28,212 subjects passing the iterative QC procedure were smoothed with a 4mm full-with-half-maximum gaussian kernel in SPM12⁷⁸, concatenated across all subjects and multiplied with a binary mask to isolate voxels located in the cerebellar cortex. This cerebellar cortical mask was constructed by multiplying a binary mask containing all 28 cerebellar lobules of the SUIT-atlas with the thresholded (at a value of 0.1) mean (unsmoothed) grey matter segmentation across all 28,212 participants.

The smoothed, concatenated and masked grey matter maps were then subjected to orthogonal projective non-negative matrix-factorization (OPNMF)²⁶, in order to derive data-driven parcellations of regional cerebellar grey matter volume.

Non-negative matrix factorization (NNMF) is a blind source separation technique that allows structural brain networks to be described in a hypothesis-free, data-driven way by identifying patterns of covariation in the data. In contrast to alternative techniques, such as principal component analysis and independent component analysis, which yield components with both positive and negative weights that are often difficult to interpret, NNMF produces a sparse, positive-only, parts-based representation of the data. Importantly, NNMF has previously been proven effective in estimating covariance patterns in neuroimaging data while providing an easier interpretation of the results than other matrix decomposition techniques such as principal component analysis (PCA) or independent component analysis (ICA)^26–30.

Briefly, NNMF decomposes an input matrix (voxels×subjects) into two matrices; a component matrix W (voxels×k) and a weight matrix H (k×subjects) where k is the number of components that needs to be specified by the user. Here, we used an implementation of orthogolnal projective non-negative matrix factorization previously used in a number of publications^26–30 and downloaded from https://github.com/asotiras/brainparts. Given the large number of participants, we used the opnmf_mem.m function, which has been optimized for high-dimensional data. The function was run with default parameter settings, exept for maximal number of iterations, which was increased to 200k, in order to ensure full convergence across all tested model orders.

Model order selection

Since the resulting parcellations are highly dependent on the requested number of components (or model order) specified, we tested model orders ranging from 2 to 30 in steps of 1, as well as from 30 to 100 in steps of 10. Binarized winner-winner-takes-all (i.e. assigning each voxel to the NNMF component with the highest loading) maps of the resulting decompositions projected onto a flattened representation of the cerebellar cortex are shown in Supplementary Figure 1.

Resulting NNMF decompositions were then evaluated based on two criteria; 1) how much of the variance in the original input data a given NNMF solution (i.e., component maps and subject weights) could explain and on 2) how reproducible the component maps were. As a metric of change in explained variance we used the change of the Frobenius norm of the reconstruction error. With increasing model orders the variance explained will always increase and the reconstruction error decrease, but if the decrease in the reconstruction error (or gradient) levels off, this indicates that the intrinsic dimensionality of the data might have been approximated (and that the subsequent increase in explained variance can largely be attributed to fitting random noise in the input data). In order to assess reproducibility, we split the full dataset into two equal sets (matched with respect to scanner site, n = 14,105 and 14,107, respectively), and ran NNMF on each split-half sample. For each set of independent NNMF parcellations, we computed two reproducibility indices. First, for each model order we matched components across split-half runs using the Hungarian algorithm⁷⁹, before computing the spatial correlations between matched component maps, and extracting the median correlation across all matched components as our first reproducibility index. For our second reproducibility index, we first computed one overall – and categorical – component map using a “winner-takes-all”-approach, i.e., assigning each voxel to the NNMF component with the highest loading. Next, we calculated the adjusted Rand index (ranging from 0 to 1, with higher values indicating greater similarity⁸⁰, across the two categorical parcellation maps for each model order as our second metric of reproducibility.

As can be seen in Supplementary Figure 2A, for lower model orders, increasing the model order resulted in a sharp decrease in the reconstruction error, indicating that models with more components resulted in a significantly better fit to the data. However, after reaching model orders between 15 and 30, the incremental improvement in fit from adding another component appeared to level off. As expected, the reproducibility results showed a largely inverse pattern, with both reproducibility indices decreasing with increasing model orders (Supplementary Figure 2B). Of note, for model orders up to 8, median spatial correlations across all matched components was above .99, indicating almost identical parcellations derived from the two independent samples. However, even for a model order of 100, the median pairwise correlation was still pretty high (.85), with 60% of components showing pairwise spatial correlations above 0.8, suggesting a reasonable level of reproducibility even for the most fine-grained parcellation. Our reproducibility index for the categorical parcellations showed very similar results; the adjusted Rand index remained above .9 for model orders up to 8, and then decreased to 0.56 at a model order of 100.

Given that the change in reconstruction error appeared to stabilize between 15 and 30 components, indicating that the intrinsic dimensionality of the data had been approximated, we searched for the most reproducible parcellations within this range. Since the reliability estimates for model orders of 16 and 23 were very similar, with the 23-component solution explaining more of the variance in the original data, we used the 23-component parcellation for all further analyses.

Adjustment and rank-order normalization of anatomical indices

Prior to being subjected to genome-wide association analyses, we adjusted all anatomical indices for effects of age, sex, estimated total intracranial volume, scanner site, 40 genetic population components, genetic batch and mean Euler number (i.e, an index of MRI image quality³¹, averaged across hemispheres), using general additive models. Finally, all adjusted anatomical indices were inverse rank normalized⁸¹.

Pre-processing of genetic data

For all genetic analyses we made use of the UKB v3 imputed data, which has 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. We first selected White Europeans, as determined by a combination of self-identification as ‘White British’ and similar genetic ancestry based on genetic principal components, that had undergone the neuroimaging protocol. We then filtered out individuals with more than 10% missingness, removed SNPs with low imputation quality scores (INFO <.5), SNPs with more than 10% 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,238 SNPs. After estimating the genetic relationship matrix (GRM) using genetic complex trait analysis (GCTA³³), we finally removed 910 participants defined as close relatives using a threshold of 0.05 (approximately corresponding to 3^rd cousins).

Heritability estimation and genetic correlation analyses

SNP-based heritability estimates for all morphological features – as well as the pairwise genetic correlations between cerebellar features - were estimated using genetic complex trait analysis (GCTA³³).

Univariate genome-wide association analyses

Univariate analyses of total cerebellar grey matter volume and the ten cerebral comparison phenotypes were conducted using Plink v1.9.

Genetic correlation analyses across discovery and replication samples and with previous published GWASs

Genetic correlation analyses across different samples were conducted using LD-score regression³².

Multivariate genome-wide association analyses

For our main analysis, we used a recently developed multivariate analysis method (MOSTest³⁶), to conduct a multivariate genome-wide association (GWA) analysis on cerebellar morphological features. MOSTest identifies genetic effects across multiple phenotypes, yielding a multivariate GWAS summary statistic across all 23 features, and provides robust (permutation based) test statistics. For mathematical details of the implementation, see van der Meer et al. (2020)³⁶, for details on the software implementation see github.com/precimed/mostest. MOSTest has been extensively validated in the original methods paper, including simulations and comparisons with other methods that have confirmed its solid discovery performance as well as an order of magnitude shorter runtime compared to other tools³⁶. For comparison to standard univariate approaches, we also performed univariate GWASs (extracted from the univariate stream of MOSTest³⁶ and identical to results from analyses using Plink).

Multivariate replication analysis

To ensure that not only single locus lead SNP associations replicate, but that also the multivariate pattern of these associations are consistent in the discovery and replication sample, we implemented a multivariate replication procedure established in Loughnan et al.³⁹. In brief, for each locus lead SNP identified in the multivariate analysis in the discovery sample, this procedure derives a composite score from the mass-univariate z-statistics and tests for associations of the composite score with the genotype in the replication sample (for mathematical formulation see Loughnan et al.³⁹). 12 of the 351 locus lead SNPs could not be tested as they were not available in the replication sample after QC. For the remaining SNPs we report the percent of loci replicating at P < 0.05, the percent remaining significant after Bonferroni correction for 339 conducted tests, and the percent of lead SNPs showing the same effect direction.

Multivariate cerebral comparison phenotypes

To compare key multivariate cerebellar results with other multivariate brain phenotypes, we downloaded summary statistics from two recent studies on cerebrocortical⁴⁴ and hippocampal⁴⁵ regional morphology. These comparison summary statistics were next analyzed using FUMA³⁷ as described for the main cerebellar results below.

Locus identification and SNP annotation

To identify genetic loci we uploaded summary statistics to the FUMA platform v1.4.1³⁷ Using the 1000GPhase3 EUR as reference panel, we identified independent SNPs as SNPs below the significance threshold of P < 5e−8 that were also in linkage equilibrium with each other at r² < 0.6. For each independent SNP, all candidate variants are identified as variants with LD r2≥ 0.6 with the independent SNP. A fraction of the independent significant SNPs in approximate linkage equilibrium with each other at r² < 0.1 were considered as lead SNPs. For a given lead SNP, the borders of the genomic locus are defined as min/max positional coordinates over all corresponding candidate SNPs. Finally, loci are merged if they are separated by less than 250kb.

FUMA further annotates associated SNPs based on functional categories, Combined Annotation Dependent Depletion (CADD) scores which predicts the deleteriousness of SNPs on protein structure/function⁵⁴, RegulomeDB scores which predicts regulatory functions⁸³; and chromatin states that shows the transcription/regulation effects of chromatin states at the SNP locus⁸⁴. For all these analyses, we used default FUMA parameters.

Genome-wide gene-based association and gene-set analyses

We conducted genome-wide gene-based association and gene-set analyses using MAGMA v.1.10⁴⁹ (http://ctg.cncr.nl/software/magma). MAGMA performs multiple linear regression to map the input SNPs to 19,190 protein coding genes and estimates the significance value of that gene. Genes were considered significant if the P value was <0.05 after Bonferroni correction for 19,190 genes. The same procedure was used for MAGMA analysis of summary statistics for the three multivariate cerebral comparison phenotypes (cerebrocortical thickness and surface area⁴⁴ and regional hippocampal volumes⁴⁵).

MAGMA gene-level statistics were next used as input to gene-property and gene-set analyses in MAGMA. Gene-property analyses test for associations between tissue specific gene expression profiles and disease-gene associations. The gene-property analysis is based on the regression model: Z ∼ β₀ + E_tβ_E + Aβ_A + Bβ_B + ɛ, where Z is a gene-based Z-score converted from the gene-based P-value, B is a matrix of several technical confounders included by default (e.g., gene size, gene density, sample size), E_t is the gene expression value of a testing tissue type and A is the average expression across all tissue types in a data set (ensuring a test of expression specificity). We performed a one-sided test (β_E > 0) which is essentially testing the positive relationship between tissue specificity and genetic association of genes.

We tested associations with two regional brain gene expression datasets (Allen Human Brain Atlas³⁵ and GTEx) and one developmental brain gene expression dataset (BrainSpan). For extraction and processing of gene expression data, see below.

Extraction and processing of gene expression data

Allen Human Brain Atlas: Regional microarray expression data were obtained from 6 post-mortem brains (1 female, ages 24.0--57.0, 42.50 +/− 13.38) provided by the Allen Human Brain Atlas³⁵. Data were processed with the abagen toolbox (version 0.1.3) ³⁴ using two volumetric atlases; 1) the binarized 23-region NNMF-derived parcellation of the cerebellar cortex; and 2) a modified version of the Desikan atlas were ROIs were merged to construct 9 bilateral regions: cerebellum, cerebral cortex, pallidum, caudate, putamen, thalamus, amygdala, nucleus accumbens and hippocampus.

First, microarray probes were reannotated using data provided by Arnatkeviciute, Fulcher and Fornito⁸⁵; probes not matched to a valid Entrez ID were discarded. Next, probes were filtered based on their expression intensity relative to background noise⁸⁶, such that probes with intensity less than the background in >=50.00% of samples across donors were discarded, yielding 31,569 probes. When multiple probes indexed the expression of the same gene, we selected and used the probe with the most consistent pattern of regional variation across donors (i.e., differential stability⁸⁷), calculated with:

where p is Spearman’s rank correlation of the expression of a single probe, p, across regions in two donors B_i and B_j, and N is the total number of donors.

Here, regions correspond to the structural designations provided in the ontology from the AHBA. The MNI coordinates of tissue samples were updated to those generated via non-linear registration using the Advanced Normalization Tools (ANTs; https://github.com/chrisfilo/alleninf). Samples were assigned to brain regions in the provided atlas if their MNI coordinates were within 2 mm of a given parcel. To reduce the potential for misassignment, sample-to-region matching was constrained by hemisphere and gross structural divisions (i.e., cortex, subcortex/ brainstem, and cerebellum, such that e.g., a sample in the left cortex could only be assigned to an atlas parcel in the left cortex⁸⁵). All tissue samples not assigned to a brain region in the provided atlas were discarded. \n\nInter-subject variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function⁸⁸:

where (x) is the median and IQR_z is the normalized interquartile range of the expression of a single tissue sample across genes. Normalized expression values were then rescaled to the unit interval:

Gene expression values were then normalized across tissue samples using an identical procedure. Samples assigned to the same brain region were averaged separately for each donor and then across donors, yielding two regional expression matrices with 23 and 9 rows, corresponding to brain regions, and 15,631 and 15,633 columns, corresponding to the retained genes. Prior to inclusion in MAGMA gene property analyses, we converted gene names for the modified Desikan atlas to ENSMBL IDs, and calculated the mean expression value across tissue types (in order to include this as a covariate in MAGMA analyses testing for gene expression specificity), resulting in a 10 (regions) by (15,490) gene expression matrix.

Analysis of Human Dating Genome Data

The Atlas of Variant Age for chromosomes 1-22 was downloaded from the Human Genome Dating (HGD) website: https://human.genome.dating/. This atlas contains more than 45 million SNPs which has been assigned dates of origin based on a recombination clock and mutation clock applied to two large-scale sequencing datasets (1000 Genomes Project⁸⁹ and The Simons Genome Diversity Project⁹⁰), with no assumptions made about demographic or selective processes⁴⁶. The current study used the median joint age estimates from both clocks when analyzing SNPs present in both datasets in combination (i.e., 13,694,493 SNPs).

After merging dated SNPs with the 40,405,505 SNPs also present in the Haplotype Reference Consortium reference data (to add minor allele frequencies, MAF) and removing 715,083 (5.2%) SNPs with missing MAF values as well as the very few (14,549, 0.1%) SNPs dated older than 2 million years, 12,960,066 SNPs remained.

548 of these dated SNPs were matched to lead SNPs linked to cerebellar morphology (defined as being in mutual LD at an r² threshold < 0.1), hereafter referred to as cerebellar-SNPs. Partially because very rare variants (Minor Allele Frequency/MAF < 0.005) had been removed prior to the multivariate GWAS analysis, MAF was not equally distributed between these cerebellar-SNPs and the full range of dated SNPs in the HGD dataset. Importantly, MAF has been shown to be systematically related to the estimated age of SNPs, with a higher proportion of low-MAF SNPs in more recent than in older time-bins^46,47. Consequently, following the analysis approach established by Libedinsky et al. ⁴⁷, we first determined the MAF-distribution of cerebellar-SNPs (across eight bins: <0.0001; 0.0001-0.001; 0.001-0.01; 0.01-0.1; 0.1-0.2; 0.2-0.3; 0.3-0.4; 0.4-0.5) and selected random set of 2,892,270 SNPs from the HGD dataset that were matched to the cerebellar-SNPs in terms of MAF-bin distribution. Similar MAF-matched HDG subsets were created for lead SNPs asociated with the two cerebrocortical comparison phenotypes, i.e., regional cerebrocortical surface area and thickness (862 and 714 lead SNPs, respectively).

For statistical inference we constructed null models (separate for each brain phenotype) by randomly drawing sets of SNPs (of equal size to the number of phenotype-linked lead SNPs) from the MAF-matched HGD-datasets and computing the histograms of estimated dates from 0 to 2 million years ago (divided into 100 bins of 20.000 years) over 10,000 iterations. From these null models we extracted bin means as well as significance thresholds (defined as the upper and lower 99.95^th percentile of the null model (i.e. corresponding to a two-tailed threshold of 0.05 Bonferroni corrected across 100 bins).

For the validation analyses, we ran the same analyses based on the larger number (range across phenotypes: 1574-2883) of individual significant SNPs (defined as being in mutual LD at an r² threshold < 0.6).

Positional and eQTL mapping of SNPs to plausible causal genes

In addition to using MAGMA, we also mapped candidate SNPs to plausible causal genes using two complementary gene mapping strategies implemented ion FUMA³⁷: 1) Positional mapping of deleterious SNPs (defined as having a CADD-score > 12.37) and 2) eQTL-mapping of SNPs previously shown to alter gene expression in cerebellar tissue (from the BRAINEAC and GTEx v8 databases). These analyses were run with default FUMA parameters. For the three multivariate cerebral comparison phenotypes (i.e., cerebrocortical thickness, cerebrocortical surface area and hippocampal regional volumes), we employed identical gene mapping procedures to our cerebellar morphology results, except for the tissues chosen for eQTL mapping (cerebrocortical and hippocampal, respectively).

Gene set analyses using lists of mapped genes

All gene set analyses using mapped genes were conducted using the hypeR R-package ⁹¹. This package implements the hypergeometric test (also known as Fisher’s exact test), which assigns a p-value to gene-set overlaps given gene set sizes and the number of background genes. This R-package also contains functions for downloading and/or formatting gene sets. The following gene sets were accessed using hypeR: “Allen_Brain_Atlas_up” (regional overexpression in the Allen Muse Brain Atlas), “MGI_Mammalian_Phenotype_Level_4_2021”, “DisGeNet” (all from the enrichR platform. In addition, we downloaded sets of genes with regional overexpression in the Allen Human Brain Atlas from the Harmonizome platform⁹² and accessed a list of genes mapped to human accelerated regions⁹³. For all gene-set analyses we employed Bonferroni-correction by dividing the p-value threshold of 0.05 by the number of gene sets included in each analysis.

Genetic overlap between cerebellar morphology and brain disorders

We accessed GWAS summary statistics for attention deficit hyperactivity disorder (ADHD)⁹⁴, autism spectrum disorder (ASD)⁹⁵, bipolar disorder (BIP)⁹⁶ and major depressive disorder (MDD)⁹⁷ from the Psychiatric Genomics Consortium. In order to avoid sample overlap, for MDD we used summary statistics based on a sample with the UK Biobank participants removed. 23&me participants included in the original MDDGWAS were also excluded, since these data are not freely available). Finally, we included data from a recent study of schizophrenia (SCZ)⁹⁸. Shared variants associated with cerebellar morphology and each of the above-mentioned brain disorders were identified using conjunctional FDR statistics (FDR < 0.05)^99,100. In contrast to genetic correlation analysis, conjunctional FDR does not require effect directions and can therefore be applied to summary statistics from multivariate GWAS, which do not contain effect directions. Two genomic regions, the extended major histocompatibility complex genes region (hg19 location Chr 6: 25119106–33854733) and chromosome 8p23.1 (hg19 location Chr 8: 7242715–12483982) for all cases and APOE region for ASD, were excluded from the FDR-fitting procedures because complex correlations in regions with intricate LD can bias FDR estimation. We further controlled for spurious enrichment by calculating all conditional Q-Q plots after random pruning averaged over 500 iterations. At each iteration, one SNP in every LD block (defined by an r2 >0.1) was randomly selected and the empirical cumulative distributions were computed using the corresponding p-values. Finally, we submitted the results from conjunctional FDR to FUMA v1.3.7³⁷ to annotate the genomic loci with conjFDR value < 0.10 having an r² ≥ 0.6 with one of the independent significant lead SNPs. Genetic correlations between univariate results for the 23 cerebellar features, the 10 comparison brain phenotypes, and each of the five mental disorders were computed using LD-score regression as described above.