Methylenetetrahydrofolate Reductase Gene Variants Confer Potential Vulnerability to Autism Spectrum Disorder in a Saudi Community

Plain Language Summary

Autism spectrum disorder (ASD), a group of an early childhood-onset neurodevelopmental disorder, is characterized by the various degree of abnormal language/communication and social impairments and restricted, repetitive behaviors and interests. The current study focuses on one specific gene that the authors are interested in. Multiple genes, involved in the pathogenesis of ASD with strong genetic impact, have been identified in different cultures and ethnic peoples with ASD. The heterogeneity and the interaction effects of genetic variants could play a considerable role in such complex multifactorial disorders.

Introduction

Autism spectrum disorder (ASD) is a complex early-onset neurodevelopmental disorder characterized by struggles in social relationships, deficiency in language and speech, and stereotypical behaviors.¹ Several reports discuss clinical associations, biochemical testing, and therapeutic management of ASD cases in the Saudi population, but very few reports have been published regarding genetic variations in this population.^2–5 Due to increasing prevalence, ASD brings substantial economic and emotional burdens to affected families and societies.⁶

Genetic and environmental factors may play a substantial role in the development of ASD.^7–11 Family, and twin studies provide strong evidence supporting the contribution of genetics in the development of the disease.^7,12 Recent studies have revealed that the genes involved in the folate pathway may be risk factors for autism.^13–15 Several reports have also confirmed that mutations or polymorphic variations of some genes significantly increase susceptibility to ASD.¹⁶ Some of these genes are serine hydroxymethyltransferase (SHMT1; MIM 182144), methylenetetrahydrofolate reductase (MTHFR; MIM 607093), and methionine synthase reductase (MTRR; MIM 602568).

The MTHFR gene has been previously investigated for potential genetic associations with ASDs in Caucasians, South Asians, and North Africans.^17–20 The MTHFR catalyzes the conversion of 5,10-methylenetetrahydrofolate into 5-methyltetrahydrofolate, a substrate for the remethylation of homocysteine to methionine with the subsequent synthesis of S-adenosylmethionine. Ramaekers et al²¹ found that low levels of 5-methyltetrahydrofolate in the spinal fluid of children who had normal neurodevelopment until age four months to six months were associated with subsequent neurological regression. The addition of folinic acid as a dietary supplement corrected the symptoms. The observed favorable response to folinic acid further supports a central role for methylation in at least some developmental disorders.²¹

Polymorphisms in genes related to folate metabolism have been examined individually or in combination in the previous studies,^{17–19,22–24} but the results have been ambiguous or inconclusive. Furthermore, the polymorphisms in these genes have not been well-defined in the Saudi population. Many polymorphisms in the MTHFR gene have been identified. Among them, the 677C>T and 1298A>C loci are two damaging polymorphisms that have been confirmed to reduce enzyme activity.²⁵ It has been reported that the 677C>T SNP decreases MTHFR enzyme activity in variant homozygotes (T/T) by 70% and in heterozygotes (C/T) by 30%;^26–28 the 1298A>C SNP can also reduce this activity, though to a lesser extent.²⁵ The 677C/T-1298A/C combined genotype has been of interest because it has been associated with lower enzyme activity and higher plasma homocysteine concentrations when compared with heterozygosity for either variant.^29–31

We hypothesized that the 677C>T rs1801133 and 1298A>C rs1801131 SNPs in the MTHFR gene might be linked alone or in combination with the risk of ASD. The purpose of the present study was aiming to investigate the vulnerability of the rs1801133 and rs1801131 genetic loci to the risk of childhood ASD in a Saudi community.

Methods

Ethics Statement and Participants

Subjects have given their written informed consent, and the Institutional Biomedical Ethics Committee-Umm Al-Qura University approved the study protocol on human research based on the National Committee of Biomedical Ethics (http://bioethics.kacst.edu.sa/About.aspx?lang=en-US). The study was conducted among Saudi simplex individuals diagnosed with ASD, selected from neuropsychiatric clinics in the Western region of Saudi Arabia (including Jeddah, Mecca, and Taif), and among healthy controls without any clinical history of mental disorders, behavioral illnesses, or epilepsies.

Inclusion/Exclusion Criteria

Each diagnosis of ASD was made by a psychiatrist, neuropsychologist, or developmental pediatrician. ASD was diagnosed following the criteria of the Diagnostic and Statistical Manual of Mental Disorders (DSM-V), based on parents’ interviews, clinical observation, individual medical records, and family information.³² All cases met the minimum scores of the autism diagnostic interview-revised (ADI-R) and autism diagnostic observation schedule-generic (ADOS-G) diagnostic instruments.^33,34 We evaluated the childhood autism rating scale (CARS) to measure the severity of behavioral symptoms (a mild/moderate case with 31–36 score and a severe with 37–60 ASD).³⁵ Also, we used the Intelligence quotient (IQ) scale to assess cognitive function based on the individuals’ ages, clinical conditions, and non-verbal tests. Cases who diagnosed clinically with obsessive-compulsive disorder, attention deficit hyperactivity disorder, or any neuropsychiatric disorder were excluded. Also, we excluded cases suffering from a genetic disorder, e.g., fragile X syndrome, and microdeletion chromosomal anomalies. We selected healthy children as a control group with no family history of mental disorders, behavioral illnesses, or epilepsies‒ from the same geographical western region as the cases.

DNA Isolation

Genomic DNA samples were extracted from buccal mucosa using Oragene DNA-OGR-575 kits (DNA Genotek Inc., Ottawa-ON, Canada). Briefly, the full buccal cells were lysed in the OGR-lysis buffer at 53°C. The released DNA was obtained by precipitating with ethanol and then dissolved in elution buffer.³⁶

TaqMan Genotyping Analysis

We implemented TaqMan genotyping assays (Thermo Fisher Scientific, USA) to genotype individuals for the selected SNPs of the MTHFR gene using a 7500 Fast-Dx Real-Time Polymerase Chain Reaction (PCR) System (Applied Biosystems, Life Technologies Inc., USA). Probe assay kits were obtained from Integrated Gulf Biosystems (ABI agency, Jeddah, SA). Assay IDs for the rs1801133 and rs1801131 SNPs were C__1202883_20 and C__850486_20, respectively. All DNA samples were included in the assays. We genotyped all samples twice, and the results were 100% concordant.

Bioinformatics Analysis

We used in-silico tools to test the effect of allelic variants on their respective functional proteins. The SNPs selected for this study were missense variants (MTHFR rs1801133 “c.677C>T; p.A222V” and MTHFR rs1801131 “c.1298A>C; p.E429A”). Variant effect predictor (VEP; https://www.ensembl.org/vep) was used to predict the variant’s deleteriousness via in-silico web-tools (e.g., SIFT, PopyPhen-2, and MutationTaster). For predicting splice-site variants effect, we utilized Human Splicing Finder 3.1 (www.umd.be/HSF/), and MiRBase (http://www.mirbase.org/) was used to test the effect of variants on miRNA binding sites.

Statistical Analysis

All controls were tested for the exact Hardy-Weinberg equilibrium (HWE) using the χ² test. The P-value (< 0.05) was considered as a departure from HWE. While conducting the statistical analysis for the SNPs, we considered the codominant, dominant, recessive, overdominant, and additive models of inheritance using SNPStats software (https://www.snpstats.net). Logistic regressions for genotypic distributions and allelic frequencies for ASD cases and controls were measured in terms of odds ratio (OR) and 95% confidence interval (CI). The less Akaike information criterion (AIC) value that corresponded to the minimal expected entropy was adopted to assess the best model of inheritance. The student’s t-test and the chi-square test were used to evaluate demographic and clinical characteristics, including age, gender, IQ, and CARS score using MedCalc statistical software version 16.4.3 (MedCalc Software; https://www.medcalc.org). We used the G*Power software (http://www.psycho.uni-duesseldorf.de/abteilungen/aap/gpower3/download-and-register/) to determine sample sizes to achieve 80% power for z-testing of two independent proportions. “Priori” sample sizes were estimated provided a probability of α = 0.05, and 80% of power sensitivity using allele frequencies among ASD cases and controls.

Results

Clinical and Epidemiological Profile

We selected a total of 112 unrelated individuals with ASD, and 104 controls from the Western region of Saudi Arabia (including Mecca, Jeddah, and Taif) for the association study. Forty additional eligible patients with ASD (23 patients as their parents refused to share their clinical data, plus 17 patients as incomplete their clinical data records). As cases and controls showed nearly high percentages of male gender (85.7% and 84.6%, respectively) with no significance (t = 1.1, 95% CI 0.5–2.3; P = 0.82) (Table 1), therefore, gender should no longer influence the study design as confounder. Similarly, the distribution for ASD cases and controls showed no significance for age factor (t = 0.9, 95% CI −1.3–0.5; P = 0.34). Consequently, matching on age and gender factors can be clearly shown for convenience.

Table 1 Epidemiologic and Clinical Characteristics in ASD Cases and Controls

There was a highly significant difference in IQ between cases and controls from the Western region of Saudi Arabia (t = 5.0; 95% CI, 2.9–7.1; P > 0.0001). Cases with CARS scores ≥37 were more frequent than cases with CARS scores <37 (67.9% versus 32.1%), and this difference was highly significant (χ² = 28.6; 95% CI, 22.9–47.0; P > 0.0001) (Table 1).

Associations with SNP Variants

The 677C>T rs1801133 and 1298A>C rs1801131 SNPs showed no deviation from the expected values for HWE in controls (χ² = 1.1; P = 0.29 and χ² = 0.04, P = 0.85, respectively). The ORs of allelic variants were 5.2 (95% CI, 3.1–8.9; P < 0.0001) for rs1801133 and 22.2 (95% CI, 7.9–62.3; P < 0.0001) for rs1801131 (Table 2). The T variant allele of rs1801133 was 3.6 times more frequent in cases than in controls (36% versus 10%, respectively), and the C variant allele of rs1801131 was 15 times more frequent in cases than controls (30% versus 2%, respectively) (Table 2).

Table 2 Genotype Distributions and Allele Frequencies of the 677C>T Rs1801133 and 1298A>C Rs1801131 SNPs in ASD Cases and Controls (Crude and Adjusted by Age)

The genotypic distribution of the MTHFR rs1801133 SNP was significantly different between cases and controls under all models of inheritance (P < 0.05), except the overdominant model, with the best interactive model being the log-additive model (OR = 0.19, 95% CI, 0.06–0.60; P = 0.0012). Except for in the overdominant model, all genotypes of the MTHFR rs1801131 SNP were significantly different between cases than controls, with the best comparison in the log-additive model (OR = 0.14, 95% CI, 0.02–0.83, P = 0.0008). In the codominant model, the genotypes carrying the T variant allele of the rs1801133 SNP or the C variant allele of the rs1801131 SNP were more frequent in cases than controls (14.3% versus 0.0% for T/T, 42.9% versus 19.2% for C/T, 25% versus 0.0% for C/C, and 10.7% versus 3.8% for A/C). The homozygous variants T/T of rs1801133 and C/C of rs1801131 were absent in controls in the non-additive genotypic models. To avoid confounding in genotyping analyses, we compared the genotype distributions in cases and controls in crude or age-adjusted by matching their OR’s. Table 2 showed that there were no significant differences between the crude and age-adjusted conditions.

Associations with CARS Scores

We determined the influence of SNPs on CARS scores as a quantitative phenotype using linear regression analysis (https://www.snpstats.net). Both the rs1801133 and rs1801131 SNPs in the MTHFR gene were statistically associated with CARS scores (OR = 5.1, 95% CI, 2.8–9.3; P < 0.0001, and OR = 13.2, 95% CI, 6.7–26.0; P < 0.0001, respectively) (Table 3). Genotype distributions of these SNPs were best expressed in the codominant models for CARS score <37 and >37 (OR = 0.61, 95% CI 0.08–4.59; P = 0.004, and OR = 0.03, 95% CI 0.00–0.32; P = 0.0005, respectively). Consequently, the T variant of the rs1801133 SNP and the C variant of the 1801131 SNP were more frequent in cases with CARS scores <37 than those with CARS scores ≥37. Hence, we found that the genotypes carrying the wild-type alleles in both the rs1801133 C>T and rs1801131 A>C SNPs were significantly more frequent in cases with CARS score ≥37. Figure 1 shows that the wild homozygotes 677C/C and 1298A/A were also more frequent in cases with CARS score ≥37 (52.6% and 79.0%, respectively). The homozygous 677T/T genotype was missing in cases with CARS scores ≥37, and the heterozygous 1298A/C genotype was absent in cases with CARS score <37 (Figure 1).

Table 3 Genotype Distributions and Allele Frequencies of the 677C>T Rs1801133 and 1298A>C Rs1801131 SNPs in ASD Patients According to CARS Scores (Adjusted by Age)

Figure 1 Genotypic distributions of the 677C>T rs1801133 and 1298A>C rs1801131 SNPs in cases with CARS scores <37 and ≥37.

Combined Genotypes of 677C>T and 1298A>C SNPs

As shown in Table 4, we found six combined genotypes of the 677C>T and 1298A>C SNPs in Saudi cases: C/C-A/A (n = 48, 42.9%), C/T-A/A (n = 20, 17.9%), C/T-C/C (n = 16, 14.5%), C/T-A/C (n = 12, 10.9%), T/T-C/C (n = 12, 10.9%), and T/T-A/A (n = 4, 3.6%). We found three combined genotypes of 677C>T-1298A>C in control individuals: C/C-A/A (n = 80, 76.9%), C/T-A/A (n = 20, 19.2%), and C/C-A/C (n = 4, 3.9%). The C/C-A/C combined genotype was absent in cases with ASD, but it was rare in control Saudi individuals (n = 4, 3.9%). We also conducted an additive genotype comparison of 677C/C with 1298(A/C+C/C), 1298A/A with 677(C/T+T/T), and 677(C/T+T/T) with 1298(A/C+C/C). Table 4 shows that two of the combined genotypes were significantly more frequent in cases than controls: 1298A/A with 677(C/T+T/T) and 677(C/T+T/T) with 1298(A/C+C/C) (P = 0.049 and P = 0.0006, respectively).

Table 4 Co-Segregation of Combined Genotypes of the MTHFR 677C>T Rs1801133 and MTHFR 1298A>C Rs1801131

Haplotype Analysis and Linkage Disequilibrium

The results of the case-control haplotype analysis and comparisons of individual haplotypes between groups are presented in Table 5. Among the four possible haplotypes of the rs1801133 (C>T) and rs1801131 (A>C) loci, the C-C haplotype was found to have an overall frequency of 4.79%, with a strongly significant difference between cases and controls (P < 0.0001). The case-control haplotype analysis also showed that there was a significant difference in haplotype distribution between cases and controls (P = 0.00057). The physical distance that separates the two SNPs (rs1801133 and rs1801131) within the MTHFR gene seems to be small (1902 base pairs) (https://genome.ucsc.edu). The two SNPs were found to be in relatively strong linkage disequilibrium (D` = 0.63, r² = 0.260) (Table 6).

Table 5 Haplotype Association of the MTHFR 677C>T Rs1801133 and MTHFR 1298A>C Rs1801131 SNPs with Response to Study Groups (Adjusted by Age)

Table 6 Linkage Disequilibrium Analysis of MTHFR 677C>T Rs1801133 and MTHFR 1298A>C Rs1801131 SNPs

In-Silico Analyses of SNPs

The SIFT algorithm, Align GVGD, and MutationTaster functional prediction tools suggested that the non-synonymous c.677C>T SNP (p.A222V) may be highly deleterious to protein function, with PolyPhen-2 scores suggesting that it is probably damaging (Table S1). Similarly, the prediction tools suggested deleterious effects of the c.1298A>C (p.E429A) on ASD cases, with PolyPhen2 suggesting possible damaging effects. Splicing tools found no effect of the 677C>T SNP on the function of the MTHFR gene but found the creation of an exonic ESS site with the loss of a 5 ESE site (i.e., new acceptor site) for 1298A>C.

Discussion

This case-control study was the first investigation of associations between SNPs in a related-folate metabolism gene and ASD cases in a Saudi community. Our data showed clear associations between the 677C>T rs1801133 and 1298A>C rs1801131 SNPs within the MTHFR gene and ASD risk. These associations were statistically significant under all genetic models in this study.

Results of earlier studies in other populations have been inconclusive regarding the association between the rs1801133 and rs1801131 SNPs and ASD risk. Reports have provided no evidence of an association between the rs1801133 SNP and autism or ASD risk in diverse Egyptian,²⁰ Brazilian,³⁷ Turkish,³⁸ Romanian,³⁹ and Korean⁴⁰ populations (Table 7). In our Saudi cases, the significant association between the rs1801133T variant allele and ASD is consistent with previous results in American,^18,41,42 Canadian,¹⁷ Chinese Han,¹⁹ Indian,⁴³ and Egyptian⁴⁴ populations (Table 7).

Table 7 Correlation Between SNP Genotypes and Allelic Frequencies of MTHFR 677C>T and ASD Among Different Ethnic Populations

Of six studies in the literature, four showed significant associations between the allelic variants and genotypic distributions of the 1298A>C rs1801131 SNP and ASD.^20,40,41,44 The other two studies showed no associations between this SNP and ASD^17,18 (Table 8). James et al⁴² reported that several genes were linked to folate/homocysteine pathways in ASD cases but did not report any significant associations with MTHFR 677C>T and MTHFR 1298A>C when these SNPs were assessed separately. When these two SNPs were analyzed together, however, they were observed to contribute to increased risk for ASD.³⁹ Like James et al⁴² Paşca et al³⁹ could not find any association in a Romanian study (39 cases/80 controls) between the 677C>T variant allele and ASD. However, a significant association was found in an Indian study (138 cases) between the variant allele of 677C>T and ASD.¹⁸ Despite having strong significant associations with ASD in some populations, the rs1801133, and rs1801131 SNPs showed deviation from HWE in ASD cases and controls in those same populations.⁴¹ In a small population, the sampling of gametes and fertilization to create zygotes can cause a random error in allele frequencies, which results in a deviation from HWE (with a more significant deviation in small sample sizes and smaller deviation in large sample sizes). However, small sample sizes or inconsistencies in HWE within a study design might lead to possible biased results.

Table 8 Correlation Between SNP Genotypes and Allelic Frequencies of MTHFR 1298A>C and ASD Among Different Ethnic Populations

Furthermore, admixture among different geographical populations might increase genetic variations and perhaps create new genotypic combinations within non-isolated (or non-native) populations.⁴⁵ Thus, genetic variations among some Middle Eastern individuals (e.g., Barbarians in North Africa, Kurdish, Upper Egyptian)^46–49 should be handled with caution, as increased consanguinity, reproductive isolation, and admixture with native source populations (e.g., Black Africans, South Eastern Asians, Caucasians) have considerable roles in gene flow and founder effects in these populations. For instance, the frequency of the 677T allele is often reported to be high in Europeans and North Americans, to be low in East Asians and Africans, and to show geographical gradients in areas of Europe, North America, and India.^47,50–52 The 1298C allele has been found with a high frequency in East Asia, and with lesser frequencies in Europe, Africa, and North America.⁵³ Because these studies might include different ethnicities within a specific sample group, we may not be able to explain these conflicting results by ethnicity alone. However, the inconclusive results among different populations are likely due to differences in sample sizes, mean ages of cases, gender ratios for cases and controls, and sources of controls (i.e., hospital-based or population-based). Moreover, analytical, statistical methods can positively or negatively give rise to unfair results.

Our study revealed that the C/T-C/C, T/T-C/C, C/T-A/C, and T/T-A/A combined genotypes in the 677C>T and 1298A>C SNPs were absent in the control group among our Saudi population, which is consistent with results in other ethnic populations. In 13,473 Chinese women, Fan et al⁵⁴ reported that the C/T-C/C, T/T-C/C, and T/T-C/C combined genotypes were missing. The absence or low prevalence of these combined genotypes has been explained to be due to the small physical distance separating the 677C>T and 1298A>C SNPs, which may reduce the probability of a recombination event.^55,56 From a public health perspective, genetic testing, and prevention strategies (e.g., B vitamin fortification) based on genetic data for a population are needed to reduce the risk of diseases closely related to mutant genotype combinations (e.g., congenital anomalies).

ASD studies have shown conflicting results regarding the susceptibility of candidate genes to CARS scores.⁵⁷ El-Baz et al⁴⁴ reported that the heterozygotes 677C/T and 1298A/C were found equally (46.2%) in ASD cases with severe CARS scores, and no significant differences were found between severe and mild/moderate cases, according to CARS scores. However, in contrast to previous studies in Egyptian populations,^44,58 our study in the Saudi population showed a significant difference in the frequency of the 677C>T and 1298A>C SNPs when comparing CARS scores <37 and CARS scores ≥37 (P = 0.001, and P = 0.0014) in the non-additive models.

Based on the observed results of allelic frequencies in both cases and controls, our “post hoc” analysis for the selected SNPs—rs1801133 and rs1801131—revealed powers of 99.67% and 99.99%, respectively, among our 216 participants. According to our “priori” estimations, we would only need sampling sizes of 40 and 26 in both cases and controls for rs1801133 and rs1801131, respectively, for a power detection of 80%.

Pinning down the MTHFR gene polymorphisms for ASD has been difficult because of poor replication of studies. Firstly, some studies have had different admixed populations, or very small cohorts, which could have substantially decreased the strength of the overall results. Secondly, various studies have included different SNPs that were inconsistent with HWE in either controls or cases, which might have led to biased results of positive or negative associations. Thirdly, few studies in the literature have dealt with the combined genotypes or haplotypes of the 677C>T and 1298A>C SNPs. Fourthly, we initially wanted to test folate, homocysteine, and vitamin B12 levels and their association with ASD, but we could not track the majority of ASD cases for this assay and so excluded it from our study. Interestingly, a potential strength of this study would be achieved via combined genotype and haplotype frequencies in a Saudi population, including D` and r² statistics, as reference data, if validated and replicated in larger cohorts before clinical implication can be approved. Moreover, our present study applied SNPStats software for analyzing SNP genotyping under all models of inheritance, and the interactions of these genotypes with clinical parameters.

Conclusion

Our study showed clear evidence of associations between deleterious MTHFR 677C>T rs1801133 and MTHFR 1298A>C rs1801131 SNPs and ASD in a Saudi community. Under all models of inheritance, we found susceptibility in cases as compared with controls for these SNPs. In the analysis of clinical features of ASD, the 677T and 1298C variants were more frequent in cases with CARS scores <37 (i.e., mild/moderate phenotype) than cases with CARS scores ≥37 (i.e., severe phenotype). Regarding combined genotypes of 677C>T and 1298A>C, the presence of 677(C/T+T/T)‒1298A/A or 677(C/T+T/T)‒1298(A/C+C/C) was significantly associated with ASD risk. The C/C-A/C combination was absent in cases with ASD, but it was rare in Saudi control individuals. The C/T-C/C, T/T-C/C, C/T-A/C, and T/T-A/A combined genotypes were absent in our ethnic Saudi population. The global haplotype analysis showed a significant difference in haplotype distribution between cases and controls of 677C>T and 1298A>C SNP (P = 0.00057). Due to the small physical distance between the two SNPs, they were found to be in relatively strong linkage disequilibrium. Our study provides essential reference data for future genetic association studies and will be useful for government and health experts developing regional health management programs. Ongoing analyses of whole-exome sequences for Saudi patients with ASD are being carried out to identify more candidate genes and SNPs.