A common variant in myosin-18B contributes to mathematical abilities in children with dyslexia and intraparietal sulcus variability in adults.

The ability to perform mathematical tasks is required in everyday life. Although heritability estimates suggest a genetic contribution, no previous study has conclusively identified a genetic risk variant for mathematical performance. Research has shown that the prevalence of mathematical disabilities is increased in children with dyslexia. We therefore correlated genome-wide data of 200 German children with spelling disability, with available quantitative data on mathematic ability. Replication of the top findings in additional dyslexia samples revealed that rs133885 was a genome-wide significant marker for mathematical abilities (P(comb) = 7.71 × 10(-10), n = 699), with an effect size of 4.87%. This association was also found in a sample from the general population (P = 0.048, n = 1080), albeit with a lower effect size. The identified variant encodes an amino-acid substitution in MYO18B, a protein with as yet unknown functions in the brain. As areas of the parietal cortex, in particular the intraparietal sulcus (IPS), are involved in numerical processing in humans, we investigated whether rs133885 was associated with IPS morphology using structural magnetic resonance imaging data from 79 neuropsychiatrically healthy adults. Carriers of the MYO18B risk-genotype displayed a significantly lower depth of the right IPS. This validates the identified association between rs133885 and mathematical disability at the level of a specific intermediate phenotype.

Introduction

Mathematical ability is a quantitative phenotype for which substantial heritability estimates have been reported.[1, 2, 3] Only one genome-wide association study (GWAS) of mathematical performance has been performed to date.[4] This study involved individuals from the general population but failed to detect any genome-wide significant association. This may have been due to the contribution of multiple genes with low effect sizes, epistasis, gene–environment interactions or technical aspects of the applied pooling method.[5] Research has shown that among children with dyslexia, which is a specific developmental disorder of learning to read or to spell,[6] the prevalence of deficits in mathematical abilities is higher than in the general population.[7, 8] We hypothesized that if distinct cognitive processes underlying mathematical disability are preferentially affected in individuals with dyslexia, investigation of this population might limit the genetic heterogeneity, and thus increase statistical power for the identification of risk genes.

Although differing brain areas have been implicated in the performance of mathematic-related tasks across studies,[9] several lines of evidence suggest that the intraparietal sulcus (IPS) has a pivotal role in the representation of numerosity.[10, 11] The IPS, in particular its right horizontal segment, is robustly activated during the active processing of numbers,[12, 13] and studies of healthy individuals have revealed that this is independent of their exact notation.[12, 14, 15, 16, 17, 18] Specific IPS responses to non-symbolic numerical processing have been demonstrated in 4-year-old children with remarkable preservation of this response pattern into adulthood.[19] Abnormalities in IPS structure or function have also been detected in other disorders affecting mathematical skills, such as dyscalculia.[20, 21, 22, 23] Structural magnetic resonance imaging (MRI) has pointed out lower gray matter volume in the right IPS and bilateral cingulofrontal areas in 9-year olds with developmental dyscalculia[22] and lower gray matter volume in more widespread posterior brain regions in a similar age group.[24] Increased gray matter density was detected in the left IPS in 15-year-old dyscalculics with preceding preterm birth.[23] Turner syndrome in adults, which is associated with specific deficits in numerical processing and visuo-spatial functions, exhibited lower depth of the right IPS with trend length and volumetric deficits.[25]

The aim of this study was to identify genetic variants contributing to mathematical ability (addition, multiplication and numerosity judgements (NJs)). We applied a quantitative trait design, and used genome-wide data of German children with dyslexia. We then followed-up the most significant findings in different independent samples. Finally, a structural MRI study was performed to validate our genetic findings at the level of an intermediate phenotype.

Materials and methods

The study was approved by the respective ethics committees, and written informed consent was obtained from all probands or their legal representatives. All association tests were performed using quantitative arithmetic data. Numbers presented in this section refer to sample sizes before statistical quality control. The behavioral test batteries (Supplementary Methods) and sample characteristics (Supplementary Tables 1, 2) are presented online.

Individuals with dyslexia

German sample

The German sample comprised 386 children who were subdivided into: (i) initial GWAS (n=200, discovery sample); and (ii) first replication sample (n=186, replication-1). The children were recruited at the Universities of Marburg and Würzburg in Germany. A diagnosis of dyslexia was assigned when the result of an age appropriate spelling test was at least 1 s.d. lower than expected on the basis of intelligence quotient measurements.[26] In each individual, two mathematical abilities were tested (Supplementary Methods), among others.[27] Briefly, arithmetic abilities were assessed using addition and multiplication tasks. The two scores were then combined, and the mean was calculated as a measure of mathematical calculation phenotype (mathematical calculation, MC). A second task addressed NJ. For NJ, subjects were asked to make judgements on the number of objects or dots presented. In a previous study, a principal component analysis revealed that MC and NJ contributed to a common ‘basic mathematical ability factor' (BMAF).[27] This factor was recalculated for the present data, and was used as a third quantitative phenotype.

German–Austrian replication-2 sample

A total of 206 German (Munich) and 118 Austrian (Salzburg) children were recruited within the context of the NeuroDys consortium. A diagnosis of dyslexia was assigned if reading performance was 0.75 s.d. below that expected on the basis of age and intelligence quotient. Mathematical ability was assessed by MC and NJ. For the genetic association analysis, MC and NJ were Z-standardized, and the mean of both phenotypes was used as a combined measure (BMAF).

Non-affected samples

German–Austrian control sample

The ascertainment process and test batteries used for the 214 controls from Munich and the 207 controls from Austria were analogous to those used for the German–Austrian replication-2 sample. To be included as a control, the child was required to have a reading performance of 0.85 s.d. above the age norm.

The mathematical test batteries (addition, multiplication and NJ) were analogous to those used for the German–Austrian replication-2 sample.

Twin Early Development Study (TEDS) sample

This sample constituted a population-based sample. Children were recruited as previously described,[28] and each performed tests of reading, spelling and mathematical ability at the ages of 7, 9 and 10 years. The mathematics-related phenotypes were assessed by: (i) using and applying mathematics (including addition and subtraction tasks); (ii) number tasks (including NJ and counting); and (iii) the perception of shapes, space and measures. All data were ascertained according to an ordinal scale. To generate a quantitative phenotype, standardized scores were summed and re-standardized. Given the genetic continuity of quantitative traits,[2] an average of this score across the three age groups was used. In the present analyses, 1081 unrelated individuals with measurements at all three time points were randomly chosen (one from each twin pair). This sample comprised 356 individuals who had been included in the previous study by Docherty et al.[4] However, as rs133885 was not one of the individually genotyped markers used in that study, this sample can be considered an independent replication sample.

Genotyping

The genome-wide data were retrieved from a previous GWAS on dyslexia that had been conducted using the Illumina HumanHap300k BeadChips (Illumina, San Diego, CA, USA). In the replication studies, single-nucleotide polymorphisms (SNPs) for the single-marker and haplotype analyses were genotyped using the Sequenom MassARRAY iPlex Gold system (Sequenom, San Diego, CA,USA) and the SpectroTYPER software v3.1 (Sequenom).

Statistical analysis

For the GWAS data, analysis of all autosomal SNPs was performed using analysis of variance, the WG-Permer software (http://www.wg-permer.org/), and R,[29] as previously described.[30] Genome-wide significance was defined according to Dudbridge and Gusnanto (P<7.2 × 10−08).[31] The quality parameters used in the replication samples and the GWAS sample were identical. Quantitative association analyses of the replication and combined samples were also conducted with R (using the glm option in the R-library stats). The covariates were: (i) sample origin (replication-2); and (ii) sample origin, spelling and word reading (combined samples).

Haplotypic effects were first tested in the genome-wide data, using a sliding window approach for up to five SNPs (including rs133885). Nominal significant haplotypes were re-analyzed in the replication-1 sample. Haplotypic P-values and linkage disequilibrium measures were calculated using PLINK.[32] Power analysis was performed post hoc, using the Genetics Power Calculator[33] and the settings ‘no dominance', ‘singletons with no parental genotypes' and a genotypic frequency of 0.317 (corresponding to the GG genotype in CEU HapMap-individuals). Effect sizes were estimated using the R2-goodness-of-fit measure and two settings: (i) the R2 for each subsample; and (ii) the weighted average of R2 across both dyslexia- and non-affected samples, respectively (Supplementary Table 3).

Sequencing analysis

Sequencing of the first four exons of MYO18Blong was performed exon-wise, in 35 children (replication-1 sample), according to standard procedures. Briefly, PCR products were purified using Agencourt AMPure (Beckman Coulter, Krefeld, Germany) and sequenced in both directions. Sequencing products were purified by Agencourt CleanSEQ (Beckman Coulter, Krefeld, Germany) and run on an ABI Capillary Sequencer (3130XL, Applied Biosystems, Darmstadt, Germany). Data were analyzed using ChromasLite 2.0 (http://technelysium.com.au/). Detected variants were followed-up in the entire discovery sample. Primer sequences and protocols are available on request.

Expression analysis

Expression analysis for MYO18B was performed using N-terminal primers for MYO18B, and primers for the core region of MYO18B (overlapping in MYO18B and MYO18B). Primers for glyceraldehyde 3-phosphate dehydrogenase were used as housekeeping controls. Complementary DNAs were obtained from customary multiple tissue complementary DNA panels (Human MTC-panel Human I (Lot-Nr.7080213); and Fetal (Lot-Nr. 7060067, all Clontech/Takara Bioscience, Saint-Germain-en-Laye, France)). PCR products were assessed on a 1% agarose gel. MYO18B identity was confirmed by sequencing. Primer sequences and protocols are available on request.

Structural MRI analysis

A structural MRI study was performed to measure sulcal depth in the IPS. This involved automated sulcus allocation and a feature extraction pipeline. Voxel-based morphometry (VBM) was used to explore genotype-dependent differences in regional gray matter volume. The genotypic and imaging data were retrieved from an independent genetic study[34] of 100 neuropsychiatrically healthy adult subjects. All details of subject characteristics, MRI sequence parameters, image postprocessing steps and the software used for sulcus morphology analysis and VBM are provided online (Supplementary Methods, Supplementary Table 4). Briefly, sulcus morphology analysis was directed to the bilateral IPS that contains the horizontal segment (subsequently referred to as IPSmain). The bilateral central sulcus was chosen as a control structure. The average sulcus depth and the total sulcus volume of all classified sulci in each hemisphere were calculated for later use as covariates. Statistical analysis was performed using analysis of covariance, and the covariates age, gender, and the average sulcus depth or the total sulcus volume. Post hoc, the sulcus measures area surface and length were analyzed for the right IPSmain only to complete the morphometric characterization. After optimal intersubject co-registration, whole-brain VBM was performed in the same sample using the conservative cluster-based statistical inference method.

Results

Genetic studies

Initial GWAS

A Q–Q plot for the BMAF analysis is shown in Figure 1 (for MC and NJ, see Supplementary Figure 1). No significant deviation was observed. One SNP (rs133885) reached genome-wide significance with BMAF, and two SNPs (rs1399428 and rs4837521) attained this level with NJ, based on P<7.2 × 10−08 (Table 1, Figure 2, Supplementary Figure 2). None of these three markers showed significant association with dyslexia.

Figure 1

Quantile–Quantile (Q–Q) plot for the basic mathematical ability factor (BMAF). Q–Q plot of the GWAS of 200 German dyslexia patients: the distribution of expected (under the null-hypothesis) and observed χ2 values for the genotypic model are depicted. Adherence to the diagonal, which is almost perfect in the lower parts of the distribution, indicates no inflation of the statistics.

Table 1

SNPs with genome-wide significance (P<7.2 × 10–08) in the GWAS sample

Marker information				Statistical model a	P-values German sample b
SNP (phenotype)	Chr	Position (dbSNP131)	Gene (hg19)		Discovery-sample	Replication-1 sample	German sample combined
rs1399428 (NJ)	9	120 660 568	NA	genotypic	2.40 × 10-08	0.2335	7.01 × 10-06
				het./hom.	3.19 × 10-09	0.1270	1.09 × 10-06
rs4837521 (NJ)	9	120 654 461	NA	genotypic	1.53 × 10-07	0.0498	5.67 × 10-07
				het./hom.	2.03 × 10-08	0.0201c	8.46 × 10-08
rs133885 (BMAF)	22	26 159 289	MYO18B	genotypic	4.39 × 10-08	0.0045	2.75 × 10-09
				allelic	1.24 × 10-06	0.0012	9.83 × 10-09
				carrier-A	5.79 × 10-09	0.0049	8.81 × 10-10

Abbreviations: BMAF, basic mathematical ability factor; GWAS, genome-wide association study; het., heterozygous; hom., homozygous; NA, not assigned; NJ, numerosity judgement; SNP, single-nucleotide polymorphism.

Genetic models with a P-value <10−05 for the specific phenotype (discovery-sample).

P-values in bold if genome-wide significant (discovery-sample), or P<0.05 (replication-1 sample).

Risk alleles not identical between discovery and replication-1 sample.

Figure 2

Genome-wide results for the basic mathematical ability factor (BMAF). Association results of the GWAS of 200 German dyslexia patients for the genotypic model. Genomic position along the 22 autosomes is represented on the x axis; −log10(P-value) is represented on the y axis. SNPs with a nominal P-value below 10–04 are depicted as circles. The horizontal line indicates the threshold for genome-wide significance (P<7.2 × 10−08).

Replication in dyslexia samples

In the replication-1 sample, association between rs133885 and BMAF was confirmed (P=0.0049), but replication failed for the two other SNPs (Table 1). In the combined German sample, the association for rs133885 was strongest in the dominant/recessive model, with carriers of the rs133885 GG-genotype showing weaker performance in mathematical tasks than AA/AG carriers (P=8.81 × 10−10, carrier-A model, Supplementary Figure 3). Calculation of the fraction of BMAF variance that was accounted for by rs133885 yielded an estimated effect size of 9.4%. Haplotype analyses did not improve P-values for the association with BMAF, compared with rs133885 alone (Supplementary Table 5).

Genotyping rs133885 in the dyslexia replication-2 sample revealed a significant BMAF association (P=0.0446). Combining all cases resulted in a P-value of 7.71 × 10−10. At least a trend toward association was found for MC and NJ, but neither of these two measures was solely responsible for the association (Supplementary Table 6). As some of the dyslexia cases showed very low values for BMAF, an inverse rank-based transformation was performed to achieve a more normal distribution. Thereafter, the P-value remained genome-wide significant (P=8.43 × 10−10, Supplementary Table 7). The overall effect size for the association between rs133885 and BMAF in samples with dyslexia was 4.87% (Table 2).

Table 2

Association analyses for rs133885 in dyslexia, control and population-based samples

Sample	Total number of individualsa	Genotypes AA/AG		Genotype GG		P-valueb	Effect sizec
		N (Frequency)	BMAF mean±s.d.	N (Frequency)	BMAF mean±s.d.
Discovery	200	132 (0.66)	−0.31±1.55	68 (0.34)	−2.10±2.61	5.78 × 10−09	15.78%
Replication-1	183	126 (0.69)	−0.11±1.90	57 (0.31)	−0.89±1.80	4.91 × 10−03	3.63%
Replication-2	316	229 (0.72)	−0.39±0.66	87 (0.28)	−0.55±0.80	4.46 × 10−02	0.88%
All case samplesd	699	—	—	—	—	7.71 × 10−10	4.87%
German–Austrian controls	391	275 (0.70)	0.42±0.77	116 (0.30)	0.41±0.81	4.86 × 10−01	0.01%
UK TEDS	1080	741 (0.69)	0.10±0.85	339 (0.31)	0.01±0.79	4.82 × 10−02	0.26%
All control/ population samplesd	1471	—	—	—	—	7.52 × 10−02	0.007%

Abbreviations: BMAF, basic mathematical ability factor, as defined in the respective samples; N, number of observations; TEDS, Twin Early Development Study.

After quality control. No deviation from Hardy–Weinberg was observed for any of the samples. Demographic details of the samples are presented as Supplementary Material online.

One-sided, except for the German discovery and combined samples (carrier-A model). Sample origin was used as a covariate in replication-2.

Measures based on R2-goodness-of-fit coefficient. In the presence of covariates, effect size was based on the difference between the R2-coefficient in the full model, and the R2-coefficient including only the covariates in the regression model.

Analysis of combined cases and controls involved use of sample origin, spelling and word reading as covariates.

Replication independent of dyslexia

Genotyping rs133885 in the population-based TEDS sample revealed a significant P-value of 0.048 (Table 2). As the phenotypes ‘numbers' and ‘using and applying mathematics' are the most comparable to BMAF, we recalculated the association for TEDS after removing the phenotype ‘shapes, space, and measures'. The P-value remained stable (P=0.041; Supplementary Table 8). The rs133885 association in the German–Austrian control sample was in the same direction, although it did not reach significance (P=0.486, Table 2).

The SNP rs133885 maps to the coding region of the Myosin-18B gene (MYO18B, MIM607295). Two distinct MYO18B protein isoforms have been identified: MYO18Blong (2567 amino acids, SwissProt accession number Q8IUG5-1); and MYO18Bshort (2080 amino acids, Q8IUG5-2). The p.E44G substitution mediated by rs133885 is located in the N-terminus specific to MYO18Blong. Sequencing of this N-terminal region revealed the presence of three common variants, two of which (rs61734946 and rs13058434) were non-synonymous. In the discovery sample, rs61734946 was strongly associated with mathematical abilities, but this effect was driven by rs133885 (Supplementary Table 9). The SNP rs13058434 (p.P177L) showed nominally significant association with BMAF in the subgroup of children who were not-at-risk according to rs133885 genotype.

Glyceraldehyde 3-phosphate dehydrogenase expression was found in all analyzed tissues and stages. MYO18B was expressed in all tissues except brain and adult lung, and strong expression was found in skeletal muscle and heart. The core part of MYO18B was abundant in all tissues, including adult brain. However, this expression was not observed in the developing fetal brain (Supplementary Figure 4).

The results of the sulcus analysis on 79 individuals are shown in Figure 3. In carriers of the rs133885 GG-genotype, a lower average depth of the IPSmain (mean±s.d. 23.18±2.62 versus 24.70±2.05 mm; F1, 74=11.74, P=0.001) was observed after adjustment for age, gender and average right-hemispheric sulcus depth. Similar results were obtained for right IPS sulcus volume (1665±535 versus 1900±425 mm3; F1, 74=4.69, P=0.033; see Supplementary Results for surface area and sulcus length). In the left hemisphere, genotype effects were in the same direction, with a trend level for IPS depth (24.79±2.15 versus 25.41±2.14 mm; F1, 74=3.69, P=0.059); but nonsignificant for sulcus volume (P=0.167). No genotype effect was found for the depth or volume of the central sulcus in either hemisphere (P>0.4; Figure 3).

Figure 3

Morphological analysis of the intraparietal sulcus (IPS). (a) Automated assignment was performed for right IPS segmentations and central sulcus. Three representative cases are shown. Note inter-individual variance of IPS geometry. Structural measurements were made for the depth (b) and volume (c) of the hypothesized sulcus and a control sulcus (central sulcus). Analysis of covariance comparing A-allele carriers with non-A-allele carriers was performed for both measures, and was corrected for age, gender, and hemisphere averages for sulcal depth and volume, respectively. Bars show estimated means and 1 s.e.m. *P<0.05; **P<0.0063 (Bonferroni-corrected threshold).

The effect on the depth of the right IPSmain withstood correction for multiple testing (Pcorrected=0.008). Regional specificity was emphasized by an explorative analysis of all extracted sulci, which yielded no stronger association than that obtained for IPSmain (Supplementary Results). VBM analysis revealed no significant genotype effect on regional gray- or white matter volume (Supplementary Results).

Discussion

In this study, we systematically investigated genetic variants contributing to mathematical abilities. Genome-wide data were retrieved from a GWAS of dyslexia, in view of previous reports that the prevalence of mathematical deficits is higher in dyslexia compared with the general population.[7, 8] Phenotypic data on mathematical performance were available for all children, thus enabling a quantitative approach.

Although three variants were genome-wide significant, only rs133885 was successfully replicated. The other two SNPs are located on chromosome 9q33.1 (rs1399428 and rs4837521) and showed genome-wide significance with NJ. They are around 6.1 kb apart and share high LD (r=0.71, D′=0.88), indicating that they reflect the same effect. Unexpectedly, however, we failed to replicate this association in the replication-1 sample, which may be due to a number of reasons. First, children included in the GWAS were more severely affected with dyslexia compared with children in the replication-1 sample. It is therefore possible that the genetic effect of the 9q33.1 locus on mathematical ability is triggered by the severity of spelling impairment. A second possible explanation is that the two markers appeared among the top GWAS hits simply by chance, without a true contribution to the quantitative trait analyzed.

It also has to be mentioned that it is highly possible that genetic components mainly defining NJ or MC may exist. By restricting our replication study to the three top hits, we might have missed other possible causal variants. Further experiments are therefore warranted to dissect the contribution of the 9q33.1 locus and of other yet unknown loci to arithmetic abilities.

rs133885 was replicated with an estimated effect size of 9.4% in the German sample. Although this might be an overestimate because of a ‘winner's curse' effect in the GWAS data, the observation of 3.6% in the replication-1 sample independently supports a strong contribution of rs133885 in individuals with spelling disorder. The difference in effect sizes, however, may also be attributable to the different degrees of dyslexia severity in the samples. Results of a stratification analysis in the German sample support this hypothesis, as an increasing effect size was observed with increasing degree of spelling disorder (Supplementary Table 10). An even lower effect size was observed in the replication-2 sample. Although this may have resulted from simple power issues, this observation might also reflect potential differing underlying core deficits in children with spelling disorder (initial, replication-1), and those with a disorder in reading (replication-2).

Although different phenotypic measures were used for the assessment of the complex construct of mathematical abilities, a consistent pattern of findings emerged from our data: the genetic association appeared to be equally strong for MCs and NJ. Our measures of mathematical abilities (accuracy and fluency of mental addition and multiplication) are influenced by verbal as well as general educational deficits that are typically experienced by dyslexic individuals and the genetic association might reflect these particular problems. However, this is not the case for the very simple task of number comparison. It has been claimed that arithmetical abilities are built on an inherited core knowledge system that can be assessed by such simple tasks[35, 36, 37] suggesting that our findings are a direct indication for the genetic basis of this numerical core knowledge, rather than mathematical ability as such.

The aforementioned results suggest that rs133885 contributes to mathematical performance in children with dyslexia, and that a more pronounced effect is found in those with spelling disorder. On the basis of this, we investigated whether rs133885 was associated with mathematical performance independent of weak reading or spelling disabilities. Although we observed a nominally significant association in the population-based sample, no statistically significant association was observed in children with superior reading performance (unfortunately, no sample showing superior spelling performance was available). However, the observation of a trend in the hypothesized direction suggests that an effect of rs133885 might be detected in larger samples (Supplementary Table 3). We also attempted a subgroup analysis of the TEDS, involving only those individuals with a weak performance in both the mathematical and the reading–spelling tasks. However, only a few individuals met these criteria, and so this analysis could not be completed (data not shown). In total, our results give rise to the hypothesis that the association between rs133885 and mathematical abilities is enhanced in dyslexic individuals. One possible explanation therefore might be the lack of compensatory mechanisms in dyslexic individuals, which render the genetic effect more pronounced at phenotypic level and, thus, more likely to be detected in smaller samples.

The associated variant rs133885 encodes an amino-acid within MYO18Blong (p.E44G). MYO18B belongs to the family of unconventional myosins,[38] and has been implicated in carcinogenesis[39, 40] and the composition of myofibrillar structures.[41] The functionally relevant position of rs133885 suggests that the variant is causative, and this was supported by the results of the haplotype analyses. Based on recent evidence suggesting the presence of important actin-binding sites within the MYO18B N-terminus[42, 43] (Supplementary Figure 5, Supplementary Results), we hypothesized that other sequence alterations of the N-terminus might also affect mathematical performance. Although no novel variants were detected, the independent association of a second common non-synonymous variant p.P177L suggests that additional N-terminal alterations within MYO18B may contribute to mathematical abilities.

The functional consequences of rs133885 are difficult to predict at the time of writing. The expression of both MYO18B isoforms in skeletal muscle and heart is consistent with the reported function of MYO18B in myocardial structures.[41] Although substantial MYO18B-expression was observed in the adult brain, no evidence was obtained for MYO18B-expression in the developing fetal brain. However, as the fetal customary complementary DNA panel comprised tissues from 20 to 30 weeks old fetuses, the relevant developmental time point might have been missed. If MYO18B is expressed at other developmental time points, this might explain the failure to detect MYO18B expression in the fetal samples. Interestingly, one previous study demonstrated a contribution of MYO18B to cognitive phenotypes.[44, 45] In one of the largest GWAS case–control studies of schizophrenia performed to date,[45] a variant located in the first exon of MYO18B was among the most strongly associated SNPs (rs5761163). Analysis of this variant in our GWAS data revealed an association with BMAF (P=0.00546), although this effect disappeared when the analysis was conditioned on rs133885 (P=0.19).

On the functional level, the sulcus analysis of the IPS demonstrated that carriers of the associated rs133885 GG-genotype showed a lower average depth of the right IPSmain. The identification of this intermediate phenotype is well concordant with the fact that basic features of gyrogenesis and sulcal morphology are under tight genetic control.[46, 47]

The lack of differences in gray matter and white matter volume as measured by VBM is divergent from previously reported VBM comparisons between dyscalculic subjects and respective age-matched healthy controls. In 9-year-old children with developmental dyscalculia, Rotzer et al.[22] found reduced gray matter volume in the right IPS and cingulofrontal areas and Rykhlevskaia et al.,[24] in the same age class, found reduced gray matter volume in more widespread occipital, cerebellar and posterior parietal regions. Isaacs et al.[23] revealed increased gray matter density, a measure sensitive to mesoscopic gyral differences, in the left parietal lobe in dyscalculic adolescents with a history of preterm birth. Several points may underlie the divergence between these results and our result: mainly, these reports studied subjects with a clinical manifestation of dyscalculia whereas the variance of the continuous behavioral phenotype explained by the genotype in this study ranged below 10%. This suggests that also at the brain level lower genotype-dependent effects may be expected compared with alterations in subjects that have impairments classifiable as dyscalculia. Although we have directed our analysis to the IPS based on converging evidence for this area being key to the representation of numerosity, formally, our MRI analysis was explorative as no generally accepted structural or functional imaging endophenotypes of arithmetic abilities (or dyscalculia) have been defined yet. Recent GWAS imaging studies, however, point out, that even for validated imaging endophenotypes such as hippocampal volume, the strongest detectable genetic univariate effects are not necessarily stronger than typical associations with a neuropsychiatric phenotype.[48] A second reason for not detecting VBM associations may lie in the age discrepancy between our MRI sample and previous VBM reports of dyscalculia. This latter assumption is strengthened by the close analogies between our results and the report of Molko et al.[25] who performed a multimodal study on adult Turner syndrome patients with dyscalculia, finding right-emphasized IPS structural differences in the sulcal analyses but only marginal gray matter alterations. Presumably, among the structural techniques, sulcus measurements could be more sensitive to gyral folding differences than VBM, which is not a surface-based technique yet, except for Molko et al.,[25] no combined sulcus/VBM analyses have been reported to allow for in-depth cross-correlation of these techniques.

Although involvement of the IPS in numerical processing has been well established,[10] recent evidence suggests that other brain areas are also involved,[9] depending on task specificities such as the type of calculation (accurate versus approximate) and type of number representation (symbolic versus non-symbolic). To determine whether our findings were specific to the IPS, we expanded our sulcal analysis to 58 sulci from each hemisphere. In the right hemisphere, one structure (superior temporal sulcus) showed nominal significance. In the left hemisphere, six additional sulci showed nominal significance with no particular regional specificity for the IPS. Despite the lack of behavioral data to demonstrate a direct association between sulcus measures and arithmetic abilities in this same study, MRI findings support the hypothesis of a specific association between rs133885 and variability in the structure of the right IPS, with weaker left-hemispheric associations possibly related to the symbolic representation of numbers in the original arithmetic tests.[9] Although there is evidence from functional studies that non-symbolic numerical processing is represented in the IPS as early as at age 4 with continuity into adulthood,[19] a direct extrapolation of our structural results to young age is hampered by limited knowledge on the trajectories of the here studied sulcus measures over the life span in humans. Hence, although the here reported sulcal variations represent a plausible first neural correlate, further structural cortex measures and functional MRI activation studies on younger samples with genetic and phenotypic characterization are essential to understand which dimensions of arithmetic operations are impacted by rs133885.

In summary, this study identified a non-synonymous variant in MYO18B as the first genome-wide significant marker to contribute to human mathematical ability. Our data provide evidence that weak spelling ability might be a modifying factor for the strength of the observed association with rs133885. The association was also present in a sample from the general population, albeit with a smaller effect size. At the level of an intermediate brain phenotype, we investigated the association between rs133885 and structural variability of the IPS, which is a key structure in terms of numerical processing. We found that healthy carriers of the MYO18B risk-genotype exhibited reduced depth of the right IPS, thus implicating MYO18B in cognitive processes. Future studies should elucidate how genetic variation in this gene exactly relates to human mathematical performance.