Genetic biomarkers predict susceptibility to autism spectrum disorder through interactive models of inheritance in a Saudi community

Abstract

Objective: To determine whether individual or interactive single nucleotide polymorphisms (SNPs) may influence the development of autism spectrum disorder (ASD). Methods: DNA from buccal cells of 212 participants (110 cases and 102 controls) were subjected to TaqMan genotyping of the HTR2A rs7997012, HTR2C rs6318, SLC6A4 rs3813034, ANKK1 rs1800497, and BDNF rs6265 SNPs. The ASD symptoms and severity were assessed by DSM-IV criteria and CARS scores. The SNPStats software was used to determine the best interactive model of inheritance of genotypic data. Results: We found susceptibility in ASD cases when compared with controls in rs7997012 (log-additive), rs6318, and rs3813034 (overdominant) and in 1800497 and rs6265 (recessive) (P< 0.05). Heterozygosity significantly contributed to the risk of ASD for rs6318 and rs3813034 SNPs (56%, P= 0.03 and 89%, P= 0.005, respectively). The rs6318 and rs6265 SNPs were significantly associated with cases with CARS scores ≥37 (recessive) (P= 0.03 and P= 0.05, respectively). Both the rs7997012 and rs6265A variant alleles were strongly associated with ASD cases with CARS scores ≥37 (P= 0.005 and P= 0.003). Conclusions: Our study provides clear evidence of associations between all five examined biomarkers and risk for ASD. Achieving exome analyses for Saudi patients with ASD could enable to identify more genetic variants and candidate genes.

Keywords:

• autism spectrum disorder
• single nucleotide polymorphism
• TaqMan genotyping
• CARS score
• interactive model of inheritance
• Saudi community

PUBLIC INTEREST STATEMENT

Autism spectrum disorder (ASD), a group of an early childhood-onset neurodevelopmental disorder, is characterized by the various degrees of abnormal language/communication and social impairments and restricted repetitive behaviors and interests. 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.

Competing interests

The authors declare no competing interests.

1. Introduction

Several reports have affirmed that mutations or polymorphic variations of some genes significantly increase susceptibility to autism spectrum disorder (ASD) (Abrahams & Geschwind, 2008). Individuals with fragile X syndrome are usually clinically co-diagnosed with ASD (AlOlaby et al., 2017; Hogan et al., 2017; Kaufmann et al., 2017). Previous reports proposed that ASD, similarly to fragile X syndrome, may be associated with language complications. But ASD may show distinct differences based on the language profiling pattern for similar developmental levels (Thurman, McDuffie, Hagerman, Josol, & Abbeduto, 2017).

The serotonin transporter gene SLC6A4 HTT has been linked to 17q11-12, as reported in several genome-wide association studies (GWAS) for autism (Auranen et al., 2002; Strobel et al., 2007). Other serotonin-related genes, such as HTR2A, HTR3A, HTR5A, and HTR1B (encoding the serotonin 5-HT1B, 5-HT2A, 5-HT3 and 5-HTR5A receptors, respectively) (Cho, Yoo, Park, Lee, & Kim, 2007; Coutinho et al., 2007; Lappalainen et al., 1995; Orabona et al., 2009), as well as genes encoding the D1 and D3 dopamine receptors (de Krom et al., 2009; Hettinger, Liu, Schwartz, Michaelis, & Holden, 2008), have also been proposed to be associated with autism susceptibility.

Like other neuropsychiatric disorders such as schizophrenia and depression (Aas et al., 2018; Connolly et al., 2006; Corponi et al., 2018; Hashimoto et al., 2005), ASD is likely associated with the brain-derived neurotrophic factor (BDNF). The BDNF gene (MIM 113505), located at chromosome 11p14.1, shares in the activity-dependent synaptic development and survival of serotonergic neurons (Cheng et al., 2009; Daws, Munn, Valdez, Frosto-Burke, & Hensler, 2007; Nishimura et al., 2007; Ren-Patterson et al., 2005; Yangngam et al., 2014). The serotonergic system and BDNF show a role in ASD and social behaviors (Dolen, Darvishzadeh, Huang, & Malenka, 2013).

Genotyping and GWAS have shown that the neurotransmitter dopamine is associated with specific brain areas and linked with different dopaminergic functions, depending on the type of autistic disorder (Ayalew et al., 2012; Gangi, Messinger, Martin, & Cuccaro, 2016; Hettinger et al., 2012; Lee, Raygada, & Rennert, 2012). The rs1800497 single nucleotide polymorphism (SNP), located just 10 kb away from the DRD2 gene in the 3` untranslated region (UTR 3`), has been identified in exon 8 of the ANKK1 gene (MIM 608774, 11q23.2). However, this functional SNP seems to be in linkage disequilibrium with several DRD2 genetic variants (Arab & Elhawary, 2015; Doehring et al., 2009). Therefore, we propose a potential dopaminergic role in the etiopathogenesis of ASD.

Several reports in the literature discuss clinical correlates, biochemical testing, and therapeutic management in Saudi cases of ASD, but very few reports have been published regarding genetic variations in this population (El-Ansary & Al-Ayadhi, 2014; El-Ansary, Bhat, Al-Daihan, & Al Dbass, 2015; Halepoto, Al-Ayadhi, & Salam, 2014; Halepoto, Bashir, Zeina, & Al-Ayadhi, 2015; Mostafa, Bjørklund, Urbina, & Al-Ayadhi, 2016). However, several polymorphisms in the genes HTR2A (MIM 182135, 13q14.2), HTR2C (MIM 312861, Xq23), and SLC6A4 (MIM 182138, 17q11.2) have been suggested to be associated with ASD (Veenstra-VanderWeele et al., 2002; Guhathakurta et al., 2009; Hranilović et al., 2010; Kistner-Griffin et al., 2011).

We hypothesized that selected genetic variants of the HTR2A, HTR2C, SLC6A4, ANKK1, and BDNF genes were associated with the occurrence of ASD. Our study specifically explored the associations of the HTR2A rs7997012, HTR2C rs6318, SLC6A4 rs3813034, ANKK1 rs1800497, and BDNF rs6265 SNPs with risk of ASD in the Saudi community.

2. Clinical significance

Although several reports in the literature discuss clinical correlates, biochemical testing, and therapeutic management in Saudi cases with autism spectrum disorder (ASD), very few reports have been published regarding genetic variations in this population. This unreeled investigation presents clear evidence of influences of all five examined biomarkers on the clinical assessment and the severity of ASD in different patterns of inheritance. However, actual collaborative work between neonatologists, geneticists, and psychologists may promote and enhance the genetic diagnosis, and earlier intervention of ASD and facilitate the potential therapies in these cases.

3. Subjects and methods

3.1. Study population

The study was conducted among unrelated Saudi patients diagnosed with ASD and healthy controls. The case group consisted of 110 cases (96 males, 14 females) aged 5–14.5 years, selected from neuropsychiatric clinics in the Western governorates of Saudi Arabia (Figure 1). Patients with ASD were clinically diagnosed using the criteria of the Diagnostic and Statistical Manual of Mental Disorders-Fourth Edition (DSM-IV), based on parents’ interviews, clinical observation, medical records of the individuals, and family information (Battle, 2013). All cases met the DSM-IV criteria for autism and minimum scores of the autism diagnostic interview-revised (ADI-R) and autism diagnostic observation schedule-generic (ADOS-G) diagnostic instruments (Lord et al., 2000, Lord, Rutter, & Le Couteur, 1994). Moreover, the childhood autism rating scale (CARS) was utilized to measure the severity of behavioral symptoms of cases; CARS score of 30–36 for mild to moderate autism, and a score of 37–60 for severe autism (Schopler, Reichler, DeVellis, & Daly, 1980). Cognitive functions were assessed using Wechsler IQ scales, according to age, clinical condition, and non-verbal tests. Cases who clinically had obsessive-compulsive disorder, attention deficit hyperactivity disorder, or any neuropsychiatric disorders were excluded. Also, those probands with a genetic disorder (e.g. fragile X syndrome, microdeletion chromosomal abnormalities) or a positive family history of a known genetic disorder was excluded. Healthy controls (102 individuals; 88 males, 14 females) were selected as having normal range of CARS scores (unpublished data) and without a family history of mental disorders, behavioral illnesses, or epilepsies.

Figure 1. Flow chart of the number of eligible individuals, and applied methodology in the study population.

The present work was approved by the Biomedical Ethics Committee, College of Medicine, UQU (#HAPO-02-K-012) (http://bioethics.kacst.edu.sa/About.aspx?lang=en-US). Parents of all participants provided informed consent.

3.2. DNA isolation

Genomic DNA samples were obtained from buccal swaps using Oragene DNA-OGR-575 kits for non-sputter (DNA Genotek Inc., Ottawa-ON, Canada). Briefly, buccal cells were collected within 30 min and lysed with the OGR-buffer at 53°C to release the DNA by precipitation with ethanol (Elhawary et al., 2017).

3.3. TaqMan genotyping

We adopted TaqMan genotyping assays (Thermo, Applied Biosystems, USA) to genotype individuals for the selected SNPs of the HTR2A, HTR2C, SLC6A4, ANKK1, and BDNF genes (Table 1) using a 7500 Fast-Dx Real-Time PCR System (Thermo, Applied Biosystems, Life Technologies Inc., USA). To ensure the accuracy of genotyping, all DNA samples, as well as negative controls, were included in the assays. We repeatedly genotyped 10% of the samples, and the results were 100% concordant.

Table 1. Examined SNP markers

Marker name	Assay ID	NCBI SNP reference	Celera ID	SNP type
HTR2A	C_1619749_10	rs7997012	hCV1619749	Intron
HTR2C	C_2270166_10	rs6318	hCV2270166	Missense mutation
SLC6A4	C_27504840_10	rs3813034	hCV27504840	UTR 3`
ANKK1	C_7486676_10	rs1800497	hCV7486676	Missense mutation
BDNF	C_11592758_10	rs6265	hCV11592758	Missense mutation

3.4. Bioinformatics analysis

We used in silico tools to test the effect of different variants on their respective functional proteins. Among the five SNPs selected for this study, three were missense variants (HTR2C rs6318, ANKK1 1800497, and BDNF rs6265), one was in the intronic region (HTR2A rs7997012), and one was in the UTR-3` (SLC6A4 rs38113034) (Table 1). Sorting Intolerant From Tolerate (SIFT) (http://sift.jcvi.org/), PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/), MutationTaster (http://mutationtaster.org/), and Align GVGD (http://agvgd.hci.utah.edu/agvgd_input.php) were examined for predicting mutations, splicing sites, and linkage to miRNA binding sites (https://www.ensembl.org/). Human Splice Finder (http://www.umd.be/HSF3/HSF.html) and Exon Splice Finder (http://rulai.cshl.edu/cgi-bin/tools/ESE3/esefinder.cgi?Process=home) were used to examining the effect of variants on splicing. MiRBase (http://www.mirbase.org/) was used to test the impact of variants on miRNA binding sites.

3.5. Statistical analysis

Both healthy controls and cases with ASD were tested for exact Hardy–Weinberg equilibrium (HWE) using the χ² test, and a P-value <0.05 was considered a departure from HWE. We conducted the statistical analysis for the examined SNPs considering the models of inheritance—the codominant, dominant, recessive, overdominant, and additive models—using the 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 ratios (ORs) and 95% confidence intervals (CIs). The less Akaike information criterion (AIC) value that corresponded to the minimally expected entropy was adopted to assess the best model of inheritance. The t-test and the chi-square test were utilized to evaluate demographic and clinical characteristics including age, gender, IQ, and CARS score (https://www.medcalc.org). A two-tailed P ≤ 0.05 was considered statistically significant. We used the G*Power software (http://www.psycho.uni-duesseldorf.de/abteilungen/aap/gpower3/download-and-register/) to determine adequate sample sizes of two independent proportions.

4. Results

4.1. Characteristics of the study population

For the study, 110 eligible Saudi participants with ASD (14 females and 96 males; ratio of 1: 6.9) and 102 healthy controls (14 females and 88 males; ratio of 1: 6.1) were enrolled. The average age of cases was 8.05 ± 2.08 years, with no significant difference when compared with controls (t = −0.96, 95% CI, 22.9‒47.0; P= 0.34). The mean CARS score was 41.56 ± 7.20 for all cases with ASD. The cases with CARS scores ≥37 (severe phenotype) were more frequent than those having CARS scores <37 (mild to moderate phenotype) (67.9% versus 32.7%, respectively), and this difference was highly significant (χ² = 28.6; 95% CI, 22.9–47.0; P > 0.0001). Likewise, the mean IQ score was 56.1 ± 6.5 for cases and 63.4 ± 7.8 for controls, and this difference was also highly significant (t = 7.4; 95% CI, 5.4‒9.2; P < 0.0001) (Table 2). Forty additional eligible parents with ASD could not be enrolled for different reasons, and two ASD cases were excluded from the study due to their older ages (Figure 1).

Table 2. Epidemiologic and clinical characteristics in ASD cases and controls

Parameter	ASD cases N = 110	Controls N = 102	t (95% CI) (P value)^a χ^2 (95% CI) (P value)^b
Age range, years	5‒14.5	6‒17
Mean age ± SD, years	8.05 ± 2.08	7.77 ± 2.62	‒0.96 (−0.9‒0.3) (0.34)^a*
Male gender^c	96 (87.3)	88 (86.3)
Mean IQ ± SD^a	56.1 ± 6.5	63.4 ± 7.8	7.4 (5.4‒9.2) (<0.0001)^a**
CARS scores of ASD cases
CARS range  Mean CARS ± SD^a	31‒60 41.56 ± 7.20		53.3 (40.2‒43.0) (<0.0001)^a**
Mild to moderate (<37)^c  Severe (≥37)^c	34 (32.7) 76 (67.9)		39.0 (29.2‒52.8) (<0.0001)^b**
Excluded subjects	n = 2 (age >18 years) n = 23 (refused to share) n = 17 (incomplete clinical profile)

ASD: autism spectrum disorder; CARS: childhood autism rating scale; SD: standard deviation.

* No significant difference (P > 0.05).

** Highly significant difference (P < 0.001).

^a Student’s t-test. Values are mean ± SD.

^b Chi-square test to compare the CARS scores of mild-moderate cases (<37) and severe cases (≥37).

^c Number of subjects, with percentages in parentheses.

4.2. Hardy–Weinberg equilibria of SNPs

All controls were in exact HWE at the examined SNPs: HTR2A rs7997012 (P= 0.42), HTR2C rs6318 (P= 0.09), SLC6A4 rs3813034 (P= 0.71), ANKK1 rs1800497 (P= 1.00), and BDNF rs6265 (P= 1.00) (https://www.snpstats.net). Similarly, no deviations from HWE were observed in cases with ASD at three SNPs; rs7997012 (P= 1.00), rs6318 (P= 0.41), and rs6265 (P= 0.20), but not at the rs3813034 or rs1800497 SNP (P< 0.05 each).

4.3. Associations with SNP variants

Regarding allele frequencies, the ORs of the allelic variants were 0.1 (95% CI, 0.05–0.14; P < 0.0001) for rs7997012, 1.8 (95% CI, 1.2–2.7, P= 0.008) for rs6318, 1.5 (95% CI, 1.0–2.2; P= 0.047) for rs3813034, 3.6 (95% CI, 2.2–5.8, P< 0.0001) for rs1800497, and 1.7 (95% CI, 1.1–2.6; P= 0.024) for rs6265 (Table 3).

Table 3. Genotype distributions and allele frequencies of selected SNPs in ASD cases and controls (adjusted by age)

SNP	Genotype	ASD cases (N = 110)	Controls (N = 102)	Model	Comparison	Logistic regression
		n (%)^a	n (%)^a			OR (95% CI)	P value ^b	AIC ^c
HTR2A IVS2 A > G (rs7997012):
Genotype:	A/A	82 (75)	16 (15)	Codominant	A/A vs. A/G	7.0 (1.8–31.6)	< 0.0001	53.0
	A/G	28 (25)	40 (38)	Dominant	A/A vs. A/G-G/G	16.0 (4.2–64.6)	<0.0001	60.0
	G/G	0 (0)	46 (45)	Recessive	A/A-A/G vs. G/G	NA (0.0–NA)	0.008	59.4
				Overdominant	A/A-G/G vs. A/G	1.8 (0.6–6.0)	0.34	79.4
				Log-additive	—	11.1 (3.5–38.5)	<0.0001	52.5
HWE, χ^2 (P value):		1.00	0.42
Allele:	A G	192 (0.87) 28 (0.13)	72 (0.35) 132 (0.65)			1.0 0.1 (0.05–0.14)	<0.0001	NA
HTR2C 68G > C (rs6318):
Genotype:	C/C C/G G/G	40 (36) 62 (56) 8 (7)	66 (65) 24 (24) 12 (12)	Codominant Dominant Recessive Overdominant Log-additive	C/C vs. G/C C/C vs. G/C-G/G C/C-G/C vs. G/G C/C-G/G vs. G/C —	0.2 (0.07–0.8) 0.3 (0.1–0.9) 1.7 (0.3–11.3) 0.2 (0.07–0.8) 0.55 (0.2–1.3)	0.037 0.03 0.57 0.015 0.16	76.1 75.9 80.3 73.4 78.7
HWE, χ^2 (P value):		0.41	0.09
Allele:	C G	142 (0.65) 78 (0.35)	156 (0.77) 48 (0.23)			1.0 1.8 (1.2–2.7)	0.008
SLC6A4 UTR 3` A > C (rs3813034):
Genotype:	A/A C/A C/C	8 (7) 98 (89) 4 (4)	16 (16) 54 (53) 32 (31)	Codominant Dominant Recessive Overdominant Log-additive	C/C vs. A/C C/C vs. A/C-A/A C/C-A/C vs. A/A C/C-A/A vs. A/C —	0.07 (0.01–0.6) 0.08 (0.01–0.7) 2.24 (0.4–15.0) 0.14 (0.03–0.6) 0.47 (0.2–1.5)	0.007 0.005 0.37 0.003 0.16	72.8 72.8 79.6 71.9 77.6
HWE, χ^2 (P value):		<0.0001	0.71
Allele:	C A	106 (0.48) 114 (0.52)	118 (0.58) 86 (0.42)			1.0 1.5 (1.0–2.2)	0.047
ANKK1 2137C > T (rs1800497):
Genotype:	C/C T/C T/T	58 (53) 24 (22) 28 (25)	74 (73) 28 (27) 0 (0)	Codominant Dominant Recessive Overdominant Log-additive	C/C vs. T/C C/C vs. T/C-T/T C/C-T/C vs. T/T C/C-T/T vs. T/C —	0.9 (0.3–3.3) 0.4 (0.1–1.4) 0.01 (0.00–NA) 1.4 (0.4–4.8) 0.4 (0.2–0.9)	0.80 0.15 0.0015 0.62 0.02	72.5 78.6 70.6 80.4 75.2
HWE, χ^2 (P value):		0.01	1.00
Allele:	C T	140 (0.64) 80 (0.36)	176 (0.86) 28 (0.14)			1.0 3.6 (2.2–5.8)	<0.0001
BDNF 196G > A (rs6265):
Genotype:	G/G G/A A/A	60 (55) 36 (33) 14 (13)	66 (65) 32 (31) 4 (4)	Codominant Dominant Recessive Overdominant Log-additive	G/G vs. G/A G/G vs. G/A-A/A G/G-G/A vs. A/A G/G-A/A vs. G/A —	0.8 (0.3–2.7) 0.6 (0.2–1.9) 0.2 (0.02–2.2) 1.0 (0.3–3.1) 0.6 (0.3–1.4)	0.37 0.47 0.16 0.94 0.25	79.4 78.9 77.5 79.4 78.1
HWE, χ^2 (P value):		0.20	1.00
	G A	156 (0.71) 64 (0.29)	164 (0.80) 40 (0.20)			1.0 1.7 (1.1–2.6)	0.024

ASD: autism spectrum disorder; SNP: single nucleotide polymorphism; UTR 3`: untranslated region; HWE: Hardy–Weinberg Equilibrium; OR: odds ratio; CI: confidence interval.

^a Number of subjects, with percentages in parentheses.

^b P values were evaluated from logistic regression analysis after adjusting for age and gender. Bold numbers indicate statistically significant associations (P < 0.05).

^c AIC values refer to the model with the less AIC value that corresponds to the minimal expected entropy. Underlined numbers represent the best model of inheritance with the less AIC value.

In Table 4, the best interactive model of inheritance was described with the less AIC value that corresponded to the minimal expected entropy. Thus, evidence for the best interactive statistical model was log-additive for the HTR2A rs7997012 (OR = 11.1; 95% CI 3.5–38.5; P < 0.0001), overdominant for HTR2C rs6318 (OR = 0.2; 95% CI 0.07–0.8; P = 0.015), overdominant for SLC6A4 rs3813034 (OR = 0.14; 95% CI 0.03–0.6; P < 0.003), recessive for ANKK1 rs1800497 (OR = 0.01; 95% CI 0.0-NA; P = 0.0015), and recessive for BDNF rs6265 (OR = 0.2; 95% CI 0.02–2.2; P = 0.16).

Table 4. Genotype distributions and allele frequencies of selected SNPs in ASD patients according to CARS scores

SNP	Genotype	CARS score		Model	Comparison	Logistic regression
		<37 ^a	≥37 ^b			OR (95% CI)	P value ^c	AIC ^d
HTR2A IVS2 A > G (rs7997012):
Genotype:	A/A A/G	20 (56) 16 (44)	62 (84) 12 (16)	__	__	1.0 0.3 (0.04–1.5)	0.10	36.1
Allele:	A G	56 (0.78) 16 (0.22)	136 (0.92) 12 (0.08)			1.0 3.2 (1.4–7.3)	0.005
HTR2C 68G > C (rs6318):
Genotype:	C/C C/G G/G	12 (33) 16 (44) 8 (22)	28 (37) 48 (63) 0 (0)	Codominant Dominant Recessive Overdominant log-additive	C/C vs. G/C C/C vs. G/C-G/G C/C-G/C vs. G/G C/C-G/G vs. G/C —	1.8 (0.2–6.9) 0.8 (0.2–4.2) 0.0 (0.0–NA) 2.1 (0.4–10.7) 0.5 (0.1–1.9)	0.09 0.78 0.03 0.41 0.25	35.6 38.3 33.6 37.7 37.1
Allele:	C G	40 (0.56) 32 (0.44)	104 (0.68) 48 (0.32)			1.0 1.7 (0.97–3.1)	0.06
SLC6A4 3`-UTR A > C (rs3813034):
Genotype:	A/A C/A C/C	4 (11) 28 (78) 4 (11)	8 (11) 66 (89) 0 (0)	Codominant Dominant Recessive Overdominant log-additive	C/C vs. A/C C/C vs. A/C-A/A C/C-A/C vs. A/A C/C-A/A vs. A/C —	2.4 (0.1–44.5) 2.1 (0.1–38.5) 0.0 (0.00-NA) 4.9 (0.4–62.6) 0.6 (0.05–7.4)	0.27 0.61 0.13 0.21 0.68	37.7 38.1 36.1 36.8 38.2
Allele:	A C	36 (0.50) 36 (0.50)	82 (0.55) 66 (0.45)			1.0 1.2 (0.7–2.2)	0.45
ANKK1 2137C > T (rs1800497):
Genotype:	C/C T/C T/T	16 (44) 12 (33) 8 (22)	42 (57) 12 (16) 20 (27)	Codominant Dominant Recessive Overdominant log-additive	C/C vs. T/C C/C vs. T/C-T/T C/C-T/C vs. T/T C/C-T/T vs. T/C —	0.4 (0.05–2.6) 0.6 (0.1–2.8) 1.2 (0.2–8.3) 0.4 (0.06–2.5) 0.8 (0.3–2.2)	0.59 0.41 0.90 0.34 0.62	39.4 37.7 38.4 37.4 38.1
Allele:	C T	44 (0.61) 28 (0.39)	96 (0.65) 52 (0.35)			1.0 1.2 (0.7–2.1)	0.58
BDNF 196G > A (rs6265):
Genotype:	G/G G/A A/A	24 (67) 12 (33) 0 (0)	38 (51) 20 (27) 16 (22)	Codominant Dominant Recessive Overdominant log-additive	G/G vs. G/A G/G vs. G/A-A/A G/G-G/A vs. A/A G/G-A/A vs. G/A —	1.1 (0.2–6.5) 2.0 (0.4–10.6) NA (0.0–NA) 0.8 (0.1–4.3) 2.3 (0.6–8.5)	0.17 0.41 0.05 0.77 0.18	36.8 37.7 34.8 38.3 36.6
Allele:	G A	60 (0.83) 12 (0.17)	96 (0.65) 52 (0.35)			1.0 0.4 (0.2–0.7)	0.0057

ASD: autism spectrum disorder; CARS: childhood autism rating scale; SNP: single nucleotide polymorphism; UTR 3`: untranslated region; OR: odds ratio; CI: confidence interval, NA: not applicable.

^a Number of cases with CARS scores <37, with percentages in parentheses.

^b Number of cases with CARS scores ≥37, with percentages in parentheses.

^c P values were evaluated from logistic regression analysis according to CARS scores after adjusting for ages. Bold numbers indicate statistically significant associations (P < 0.05).

^d AIC values refer to the model with the less AIC value that corresponds to the minimal expected entropy. Underlined numbers represent the best model of inheritance with the less AIC value.

As for HTR2A rs7997012 A > G, five times more ASD cases than controls were carrying the A/A genotype (75% versus 15%), while no ASD cases were carrying the G/G genotype. Moreover, the heterozygous genotypes HTR2C rs6318 (G/C) and SLC6A4 rs3813034 (A/C) showed significant roles in the overdominant model when compared with homozygous genotypes (C/C-G/G and C/C-A/A, respectively). The recessive model showed that two copies of the T allele of the ANKK1 rs1800497 SNP or two copies of the G allele of the BDNF rs6265 SNP might be necessary to affect ASD risk. HTR2A A > G rs7997012 in the additive model showed that each copy of A modified the risk of ASD in cases in an additive form (2A/A + A/G to G/G).

The heterozygous genotypes of HTR2C rs6318 and SLC6A4 rs3813034 were highly increased in cases when compared with controls (56% versus 24% and 89% versus 53%, respectively). With slightly larger AIC values, the dominant models in rs6318 (OR = 0.3; 95% CI, 0.1–0.9; P = 0.03) and rs3813034 (OR = 0.08; 95% CI, 0.01–0.7; P = 0.005) suggested a role for heterozygosity in increasing ASD risk in cases. In addition, the codominant model (the most general model of inheritance) showed potential significance of the HTR2A rs7997012 A > G SNP (OR = 7.0; 95% CI 1.8–31.6; P < 0.0001), as it represented the comparison between the heterozygous genotype (A/G) and the homozygous genotype to give a non-additive risk (Table 3).

4.4. SNPs and quantitative phenotype using the CARS

We also determined the influence of SNPs on CARS scores as a quantitative phenotype using linear regression analysis (https://www.snpstats.net). Among the five examined SNPs, only HTR2C rs6318 and BDNF rs6265 SNPs were statistically associated with CARS scores <37 and ≥37 (P = 0.03 and 0.05, respectively) in the recessive model (Table 4). As for the allelic frequencies, the allelic variant rs7997012A was highly overexpressed in cases with CARS scores ≥37 compared to CARS <37 (0.92 versus 0.78; P = 0.07). The variant rs6265A allelic variant was twice expressed in cases with CARS ≥37 compared to CARS score <37 (0.35 versus 0.17, respectively). However, among all five examined SNPs, only the rs6265A allelic variant was strongly associated with cases having CARS scores ≥37 (OR = 0.4; 95% CI 0.2–0.7; P= 0.006) (Figure 2).

Figure 2. Comparison of the allelic frequencies of the examined SNPs in cases with CASRS scores <37 and ≥37. The asterisk “*” denotes the increased rs6265A variant allele in cases based on CARS scores ≥37 and <37 with a high significant difference (P = 0.006).

4.5. In-silico analyses of SNPs

The Variant Effect Predictor-genome browser (https://www.ensembl.org/) found that both the non-synonymous BDNF rs6265 SNP (c.196G > A; p.Val66Met) and HTR2C rs6318 SNP (c.68G > C; p.Cys23Ser) have had tolerated effects on SIFT scores (0.13, and 0.4, respectively) and benign effects on PolyPhen-2 scores (0.062, and 0.007, respectively). On the contrary, the Align GVGD web server found that the HTR2C rs6318 SNP activated protein function, showing the most likely and highest damage scoring (C65). The Align GVGD predicted the ANKK1 rs1800497 SNP (c.2137G > A; p.Glu713Lys) to have a profoundly damaging effect (C55) on protein function, and MiRBase found it to cause a loss of one miRNA binding site. According to the Human Splicing Finder, the HTR2A rs7997012 SNP had a proposed small effect on exon splice silencers. However, the different web servers showed no evidence of a consequence of the SLC6A4 rs3813034 SNP (UTR 3`) on protein function.

5. Discussion

This case–control study presents the first investigation of associations between ASD susceptibility and the HTR2A IVS2 A > G rs7997012 and HTR2C 68G > C rs6318 SNPs (serotonin receptor and serotonin transporter), the ANKK1 2137C > T rs1800497 SNP (DRD2 related), and the BDNF 196G > A rs6265 SNP in a Saudi community.

Overall, our results provide clear evidence of associations between the risk of ASD in the Saudi population and all five examined SNPs. Also, the higher rate of genotypic heterozygosity of the rs6318 and rs3813034 SNPs in cases suggests a possible contribution to the risk of ASD (overdominant model). Two copies of either the rs1800497T allele or the rs6265A allele might be necessary to affect the risk of ASD (recessive model). Also, each copy of A in the HTR2A rs7997012 SNP modifies the risk in an additive form, and the homozygous A/A more than doubles the risk when compared with the heterozygous A/G.

The BDNF rs6265 SNP has been reported to be associated with obsessive-compulsive disorder (Hall, Dhilla, Charalambous, Gogos, & Karayiorgou, 2003; Wang, Zhang, Zhu, Liu, & Zhou, 2015), attention deficit hyperactivity disorder (Leo et al., 2018; Xu et al., 2007; Zhang et al., 2018), anxiety-related personality traits (Lang et al., 2005; Montag, Basten, Stelzel, Fiebach, & Reuter, 2010; Wei et al., 2017), and other childhood-onset mood disorders (Kaufman et al., 2006). Despite the tolerated effect predicted by various in silico tools, the BDNF gene was previously reported to function in the hippocampal and prefrontal cortex roles involved in memory and learning of the individuals (Hariri et al., 2003). There is, however, a lack of association with BDNF levels in other neuropsychiatric conditions (Duncan, Hutchison, Carey, & Craighead, 2009; Mercader et al., 2007; Trajkovska et al., 2007; Umehara et al., 2016; Yu et al., 2008). A meta-analysis of GWAS in Sardinian people, for example, did not find associations between BDNF rs6265 and depressive disorders, although it did find associations between these disorders and both BDNF rs11030102 and NTRK3 genes (Terracciano et al., 2013). Interestingly, our study found that the variant BDNF (A) allele was three times as prevalent in ASD cases as in controls. In transmission disequilibrium testing of autistic Japanese trios, the BDNF SNP combinations were significantly associated with classical autism (Nishimura et al., 2007). However, previous studies in French, Chinese, and Korean populations failed to identify any association with rs6265 (Cheng et al., 2009; Tabagh et al., 2010; Yoo, Yang, Cho, Park, & Kim, 2014).

In contrast to our study, association studies of HTR2A in Caucasian populations have failed to detect any associations between the rs6311 (1438G > A), rs6313 (102T > C), or rs6314 (1354C > T) SNPs and autism (Veenstra-VanderWeele et al., 2002). However, studies on Korean and Indian samples have shown an association between the rs6311- rs6313- rs6314 haplotype and ASD (Cho et al., 2007; Guhathakurta et al., 2009). Another study also supported the possible involvement of these SNPs of HTR2A in the development of autism in Egyptian cases with gastrointestinal disturbances (Abdelrahman et al., 2014).

Regarding HTR2C rs6318, different models of inheritance (i.e. codominant, dominant, overdominant) showed significant associations with ASD that are inconsistent with other studies in Australian, German, French, and Brazilian populations (Hettinger et al., 2008; Paterson, 1999; Tabagh et al., 2010). However, transmission disequilibrium of HTR2C polymorphism in schizophrenic families has been reported (Paterson, 1999). Shirley et al. (2016) suggested a strong association between autism with intellectual disability and copy number variants that disrupt the gene function of HTR2C and ZEB2 (which encodes a transcription factor; MIM #605802).

The GWAS could not find any linkage between autism and the intronic region of the SLC6A4 gene but did find the linkage between autism and SLC25A12 (P < 0.00001) (Anney et al., 2010). Due to its proximity to the DRD2 gene, which is linked to comorbid neurological disorders, the missense rs1800497 SNP within the ANKK1 gene was clearly associated with risk of ASD in our study (recessive model). In contrast, Hettinger et al. (2012) showed a strong linkage between ASD (among 110 affected male patients) and the DRD2 rs1800498 C > T SNP (P = 0.0007), but not between ASD and the ANKK1 rs1800497 SNP.

ASD studies have shown conflicting results regarding associations between candidate genes and CARS scores (Halepoto et al., 2015). Meng et al. (2017) found an association serum BDNF and CARS scores, but not between the rs6265 gene polymorphism and CARS scores. Furthermore, serum BDNF was reported to be highly associated with CARS scores in a Chinese population (Zhang, Jiang, Kong, & Lu, 2014), but not in a Saudi community (Halepoto et al., 2015). Based on our results in this Saudi sample group, the rs6265 and rs6318 markers could have a significant effect on CARS scores. Also, a negative relationship was observed between 5-HT levels and the severity of the behavioral CARS score. However, the examined HTR2C rs6318 and BDNF rs6265 SNPs have potential clinical influences on the CARS assessment in the Saudi sample group.

According to our in-silico analyses of web links, the rs6265, rs6318, and rs1800497 SNPs are predicted to have seriously damaging effects on protein functions and risks for ASD phenotypes. Sometimes these SNPs appear to cause a loss of one miRNA binding site, as in the case of the rs1800497 SNP. Other times they seem to have a small effect on exon splice silencers, as in the case of the rs7997012 SNP.

Several factors may account for conflicting results in the literature regarding genetic associations with ASD. For example, several studies have included populations with ethnic admixture, different clinical criteria for ASD, or small sample sizes, which could have lessened the asset of the outcomes. Also, technologies ranging from enzymatic digestion of the PCR amplicons to high-throughput next-generation sequencing of the human genome have been used in previous studies to discriminate between genotypes within a specific SNP, or to identify many candidate genes, resulting in a broad range of false-negative or false-positive genotypic distributions. Sometimes, statistical analyses may control the significance of the outcoming data for the best interactive model of inheritance to compare single or combined genotypes among groups of interest. In the present study, we applied the SNPStats software for different models of inheritance and genetic risk for ASD.

Regarding limitations of our study, our “post hoc” analysis for the selected SNPs—rs7997012, rs6318, rs3813034, 1800497, and rs6265—revealed powers of 100%, 55.2%, 31.1%, 97.9%, and 46.5%, respectively, among our 212 participants (110 cases and 102 controls). According to our “priori” estimations, we would need sampling sizes of 12, 192, 390, 54, and 239 in both cases and controls for rs7997012, rs6318, rs3813034, rs1800497, and rs6265, respectively, for a power detection of 80%. Sometimes, overcoming limitations such as these can be difficult, but a meta-analysis may offer a feasible solution. This would enable surveying of a broad set of participants, thereby enhancing “post hoc” power and allowing a more broad-based analysis of previously available data. However, enrolling more participants within a reasonable period from a specific hospital or clinical center might have been challenging in our study. Replication of our outcomes through larger, multicenter genetic association studies will be of interest.

Suspected cases at earlier ages have been diagnosed as ASD by plasma amino acid measurements, as specific inborn errors of metabolism have been suggested to increase the risk of neuropsychological autistic behavior. (Campistol et al., 2016). Witters et al. (2016) have reported that newborn screening for cases suffering from propionic acidemia, decreased level of valine or altered valine/leucine ratio should be screened for ASD. Recently, Vargason et al. (2018) have reported that significantly increased levels of glutamate, serine, and hydroxyproline were found for the ASD cohort compared to typically developing peers (P = 0.014, 0.04, P = 0.018, respectively). However, significant partnerships between neonatologists, geneticists, and psychologists may enhance and serve the genetic diagnosis, and earlier intervention of ASD and improve potential therapies in these cases.

Of note, a Saudi study has analyzed 13 multiplex families and 27 affected individuals with ASD using Affymetrix GeneChip Mapping 250K and 6.0 arrays (Adi et al., 2015). While preparing the current manuscript, Khanzada, Butler, and Manzardo (2017) published gene analytics mapping of ASD, as well as other neuropsychiatric behavioral disorders that include impaired social interaction and communication. CHRNA7, DISC1, DRD2, FOXP2, HTR2A, MAOA, MTHFR, SLC6A3, and TPH2 genes were identified as candidate ASD genes. Analyses of next-generation sequencing using the targeted custom panel for ASD are ongoing, and some promising new data will be published in time.

6. Conclusions

Based on our results, there is clear evidence of associations between the risk of ASD in the Saudi population and all the examined SNPs: HTR2A rs7997012, HTR2C rs6318, SLC6A4 rs3813034, ANKK1 rs1800497, and BDNF rs6265. By selecting the best models of inheritance, we found susceptibility in ASD cases as compared with controls for rs7997012 (log-additive model), rs6318 and rs3813034 (overdominant model), and ANKK1 rs1800497 and BDNF rs6265 (recessive model). Heterozygosity could remarkably contribute to a risk of ASD for the rs6318 and rs3813034 SNPs.

In the analysis of clinical features of ASD, the rs6318 and rs6265 SNPs were significantly associated with cases with CARS scores ≥37 in the recessive model. Both the G allele and the A allele variants for the rs7997012 and rs6265 SNP were strongly associated with ASD cases with CARS scores. These outcomes should be handled with caution, as these SNPs do not act alone to explain such complex multifactorial disorders. Whole-exome sequencing analyses for Saudi patients with ASD are being carried out to detect more candidate genes and SNPs.

Authors’ contributions

N.A.E. and M.T.T. designed the research; N.A.E., M.T.T., I.A.S; A.D., M.R., A.A., N.Q., A.K., and E.N.E. performed the research; N.A.E., M.T.T., A.M., and E.N.E. analyzed the data; and N.A.E., M.T.T., N.B., I.A.S. and M.R. wrote the paper. Also, N.A.E. initiated the grant funding through a contract with GDRG-KACST.

Data Availability

The data sets analyzed during the current study are available from the corresponding author.

Acknowledgements

This work was funded by King Abdulaziz City for Science and Technology (KACST-GDRG-Riyadh, Saudi Arabia), under the Grant #ARP-34-44. The authors wish to thank the Faculty of Biotechnology, Modern Sciences and Arts University (MSA), 6th October City, Giza, for allowing E.N.E. to participate in silico analyses.

Additional information

Funding

This work was supported by the General Directorate for Research & Grants, King Abdulaziz City for Science and Technology, Riyadh, Saudi Arabia, under Grant #ARP-34-44. The authors wish to thank the Faculty of Biotechnology, October Modern Sciences and Arts University (MSA), Giza, for allowing E.N.E. to participate in silico analyses.

Notes on contributors

Nasser A. Elhawary

Nasser A. Elhawary (Ph.D.) is a Professor of Medical Genetics at the UQU (College of Medicine) and ASU (Faculty of Medicine-Cairo). He is conducting researches in the Molecular Genetics and Genomics laboratories (Department Head, Prof. Tayeb Mohammed T.). His research fields have focused on the molecular insights of neurogenetic disorders and neuromuscular inherited diseases as well. Extensive studies have been conducted to ASD to identify the susceptible genetic biomarkers. Our results could encourage to a broader project of exome analysis to identify more candidate genes/genetic variants. Neda Bogari (Ph.D.), Ahmad Mufti (Ph.D.), Ikhlas A. Sindi (Ph.D.), Asim Khogeer (Ph.D.), and Anas Dannoun (MSc) as a group’s principal research activity, interested in ASD studies.Mona Rashad (Prof. psych-Pediatrics), Arwa H. Arab. (Psychologist, Ph.D.), and Nermeen Qutub (Psychologist, Ph.D.) managed the diagnoses and psychological assessments. Ezzeldin N. Elhawary is an MS student (MSA, Egypt) shared the in-silico bioinformatics and molecular works.

Ezzeldin N. Elhawary

Nasser A. Elhawary (Ph.D.) is a Professor of Medical Genetics at the UQU (College of Medicine) and ASU (Faculty of Medicine-Cairo). He is conducting researches in the Molecular Genetics and Genomics laboratories (Department Head, Prof. Tayeb Mohammed T.). His research fields have focused on the molecular insights of neurogenetic disorders and neuromuscular inherited diseases as well. Extensive studies have been conducted to ASD to identify the susceptible genetic biomarkers. Our results could encourage to a broader project of exome analysis to identify more candidate genes/genetic variants. Neda Bogari (Ph.D.), Ahmad Mufti (Ph.D.), Ikhlas A. Sindi (Ph.D.), Asim Khogeer (Ph.D.), and Anas Dannoun (MSc) as a group’s principal research activity, interested in ASD studies.Mona Rashad (Prof. psych-Pediatrics), Arwa H. Arab. (Psychologist, Ph.D.), and Nermeen Qutub (Psychologist, Ph.D.) managed the diagnoses and psychological assessments. Ezzeldin N. Elhawary is an MS student (MSA, Egypt) shared the in-silico bioinformatics and molecular works.