Mining the 3'UTR of autism-implicated genes for SNPs perturbing microRNA regulation.

Autism spectrum disorder (ASD) refers to a group of childhood neurodevelopmental disorders with polygenic etiology. The expression of many genes implicated in ASD is tightly regulated by various factors including microRNAs (miRNAs), a class of noncoding RNAs ~22 nucleotides in length that function to suppress translation by pairing with 'miRNA recognition elements' (MREs) present in the 3'untranslated region (3'UTR) of target mRNAs. This emphasizes the role played by miRNAs in regulating neurogenesis, brain development and differentiation and hence any perturbations in this regulatory mechanism might affect these processes as well. Recently, single nucleotide polymorphisms (SNPs) present within 3'UTRs of mRNAs have been shown to modulate existing MREs or even create new MREs. Therefore, we hypothesized that SNPs perturbing miRNA-mediated gene regulation might lead to aberrant expression of autism-implicated genes, thus resulting in disease predisposition or pathogenesis in at least a subpopulation of ASD individuals. We developed a systematic computational pipeline that integrates data from well-established databases. By following a stringent selection criterion, we identified 9 MRE-modulating SNPs and another 12 MRE-creating SNPs in the 3'UTR of autism-implicated genes. These high-confidence candidate SNPs may play roles in ASD and hence would be valuable for further functional validation.

Introduction

Autism spectrum disorder (ASD) represents a series of neurodevelopmental disorders characterized by altered social interests, communicative deficits, and restricted and repetitive behaviors, usually with an onset before the age of three [1]. Although ASD is defined by a triad of symptoms, the levels of severity and presentation may vary among individuals, demonstrating the tremendous heterogeneity of the condition. The burden of ASD is increasing worldwide with recent prevalence estimates of 62 per 10,000 children [2]. Research in the last several years has substantiated the strong genetic basis and has further established the polygenic etiology of ASD [3,4]. For instance, a comprehensive transcriptomic analysis was performed to measure mRNA levels in post-mortem brains from autistic individuals and controls using Illumina microarrays. As a result, 444 genes were found to be differentially expressed [5]. This disease-specific, aberrant expression of genes may be caused by various factors, such as the presence of submicroscopic structural chromosomal alterations in the genome denoted as copy number variants (CNVs), or by perturbations in many regulatory mechanisms that govern gene expression [6-8]. Very recently, methylomic analysis of monozygotic twins discordant for ASD identified epigenetic variations that mediate disease susceptibility via altered gene dosage [9]. Besides gene dysfunction, studies reveal that the known etiological factors of ASD eventually converge toward the unifying theme of aberrant gene expression [5-9].

MicroRNAs (miRNAs) have emerged as key regulators for nearly 50% of protein-coding genes in the human genome, participating in all vital cellular processes [10-12]. miRNAs represent a class of short (∼22 nt) endogenous noncoding RNAs that pairs with complementary sequences called miRNA recognition elements (MREs), which are located mostly in the 3′untranslated region (3′UTR) of target mRNAs. This miRNA:MRE pairing leads to translational inhibition or mRNA destabilization by recruiting the RNA-induced silencing complex (RISC). In mammals, complete complementarity between miRNA and MRE is rare, since a minimal 6-bp match of nucleotides 2–8 at the 5′ end of miRNA (known as the seed region) is sufficient for functionality [13]. Thus a single miRNA can target approximately 200 transcripts, and more than one miRNA can act upon a single mRNA target [14]. Importantly, almost 70% of experimentally-detectable miRNAs are expressed in the human nervous system [15]. Several miRNAs orchestrate a myriad of diverse neurodevelopmental processes including neuronal cell specification and differentiation, synaptic plasticity and memory formation [16]. Accumulating evidence suggests that alterations in miRNA expression or function are associated with the cognitive deficits and neurodevelopmental abnormalities observed in autism, schizophrenia and other forms of intellectual dysfunction [17-20]. According to a recent bioinformatics study, co-expression of neural miRNAs with their target mRNAs is a common feature [21], suggesting that any interference in miRNA:MRE interactions might affect the physiological control of gene expression. Additionally, single nucleotide polymorphisms (SNPs) in the 3′UTR of genes are recognized for their ability to perturb miRNA binding by modulating existing MREs or by creating new MREs [22-25]. The pathological significance of this unique class of functional polymorphisms is gaining recognition and their biological relevance has been explored in various diseases [26,27], including neurodegenerative diseases [28,29]. To our knowledge, MRE-SNPs have not formed the basis of any previous case-control studies in autism, and their systematic assessment might extend our understanding of the complex genetic architecture of this highly heterogeneous disorder. However, the identification of ASD-associated MRE-SNPs poses a considerable challenge due to the lack of an established model that accounts for the heterogeneity of ASD and the ever-expanding modes of miRNA-mediated gene regulation.

At this point, computational approaches are still the driving force in identifying truly functional genetic variants within MREs, and numerous databases have been developed to assist in this process. Aided by a rigorous bioinformatics approach, we identified a panel of high-confidence SNPs that either modulate existing MREs or create new MREs in the 3′UTR of autism-implicated genes. We propose that alterations in the miRNA:MRE interaction caused by these MRE-SNPs could be an alternate mechanism underlying aberrant gene expression, which may contribute to predisposition, pathogenesis, inter-individual variation in gene expression, and heterogeneity of ASD (Figure 1). Future experimental validations may further clarify the effect and consequences of these candidate SNPs.

Figure 1

SNPs in MREs of autism-implicated genes may underlie ASD etiology SNPs in the 3′UTR of autism-implicated genes may either modulate existing MREs or create new MREs thereby leading to aberrant expression of these genes, which may underlie predisposition, pathogenesis or genetic heterogeneity of autism. The “DNA–brain” image is adopted from the public domain http://neurosciencenews.com/files/2013/06/dna-brain-familial-alziheimers-public.jpg, which is credited to the National Institute of Neurological Disorders and Stroke (NINDS)/National Institutes of Health (NIH) and is appropriately modified. MRE, miRNA recognition element.

Results

Study design and preliminary analyses

We followed a stepwise, integrative and stringent computational pipeline to identify and characterize MRE-SNPs in the 3′UTR of autism-implicated genes (Figure 2). A total of 484 genes were retrieved from the “Human Gene Module” of ‘SFARI Gene’, a comprehensive resource for obtaining reliable data on the genetics of ASD [30]. SNPs that modulate or create MREs across the 3′UTR of these genes were identified using MirSNP database [31]. The resultant list contained 13,945 SNPs in 468 out of the 484 autism-implicated genes. We chose MirSNP due to its high sensitivity in covering a majority of the validated miRNA-related SNPs when compared to other databases [31]. MirSNP also provides information on the plausible effect of SNPs on miRNA binding and categorizes them into one of the four functional classes. These include (i) break — SNP completely disrupts the binding of miRNA to the MRE, (ii) decrease — SNP reduces the binding efficacy of the miRNA to the MRE, (iii) enhance — SNP increases the binding affinity of the miRNA to the MRE and (iv) create — SNP creates new MREs for miRNAs. We grouped the first three functional classes as MRE modulating SNPs (MREm-SNPs) since the existing miRNA recognition sites were modified. The fourth category that created new MREs (MREc-SNPs) was considered as a separate group. Next, we checked the minor allele frequencies (MAFs) of these SNPs in Exome Variant Server (EVS) maintained by the National Heart, Lung and Blood Institute (NHLBI) Exome Sequencing Project. The EVS was chosen for several reasons. First, the database serves as a large repository for more than 10 million human SNPs identified by sequencing 15,336 genes in 6515 individuals of European American (EA) and African American (AA) ancestry [32]. Second, even SNPs with MAF as low as 0.1% can be identified owing to the larger sample size. Third, EVS also indicates the level of conservation and evolutionary constraints of the SNP loci by providing phylogenetic analysis with space/time model for conservation (phastCons) and genomic evolutionary rate profiling (GERP) scores. Finally, for many SNPs, the MAF is unavailable in NCBI’s dbSNP. Of the 13,945 SNPs we identified, it was interesting to note that only 387 spanning the 3′UTR of 181 autism-implicated genes were documented with a MAF > 0.1% in at least one population (EA or AA). This suggests that the vast majority of MRE-SNPs are either very rare or not detected in samples of the specified ancestry. Alternatively, it may indicate the action of purifying selection on MRE-SNPs [33].

Figure 2

Work flow to identify high-confidence SNPs that modulate existing MREs or create new MREs

Minimizing false positive MREs

The first step in the identification of MREm-SNPs relies on the correct prediction of functional MREs in the 3′UTR of autism-implicated genes. MirSNP database employs the miRanda miRNA target prediction algorithm [34]. Despite the high sensitivity of miRanda [35], the reliability of a single prediction method is modest. To account for this, the MREs were subjected to consistent cross-prediction by TargetScanHuman 6.2 [36,37] and RNAhybrid 2.1 [38,39]. Furthermore, we improved the stringency of the prediction by considering certain quantitative scores. As an added feature, the MirSNP database provides an miRSVR score, which was developed by supervised training on mRNA expression data from a panel of miRNA transfection experiments [40]. We considered only the top miRSVR predictions with a score of <−0.1, which have greater probability for down regulating the target mRNAs. However, MirSNP does not provide an miRSVR score for all predicted MREs. Therefore, additional care was taken to minimize the exclusion of true positives that lack an miRSVR score by considering the ‘total context score’ given by TargetScan, which effectively ranks the predicted MREs based on their local and global contextual features [37].

Obtaining high-confidence and tissue-specific datasets

The miRNA–MRE interaction involves base pairing that requires favorable hybridization energy for the formation of a functional duplex. We computed the free energy of binding of the miRNA:MRE duplex by analyzing their sequences in RNAhybrid. A maximum hybridization energy of −15 kcal/mol was used as the scoring criteria to select stable duplexes. All of these sequential filtering steps aided the identification of high-confidence functional MREs and thus, MREm-SNPs. However, the biological relevance of these SNPs depends on whether the miRNA and its target mRNAs are co-expressed in a particular tissue or cell type. Considering the neurobiological basis of ASD, we shortlisted only those MREs that served as binding sites for the brain expressed human miRNAs. This brain-specific miRNA dataset consisted of overlapping miRNAs from two independent studies that profiled the expression of miRNAs in different regions of post-mortem human brain tissue samples [41,42]. In total, 9 MREm-SNPs were identified across the 3′UTR of 17 transcripts (Table 1). The effect on miRNA binding was verified by computing the difference in hybridization energy (ΔΔG) caused by that particular SNP (see Materials and methods section). The greater the ΔΔG value is, the greater the impact of the SNP will be. The presence of the variant allele for the MREm-SNPs rs56124628 in CSNK1D encoding casein kinase 1 delta, rs28390202 in PLCB1 encoding phospholipase C beta 1, rs148637931 in PRKCB encoding protein kinase C beta, rs114517919 in SLC16A3 encoding solute carrier family 16 member 3, and rs1042589 in TPO (thyroid peroxidase) decreases the binding efficiency of their corresponding miRNAs (indicated by their lower negative ΔΔG value) to the respective mRNA transcripts and may therefore allow “leaky” translation. On the other hand, the SNPs rs112904440 in SNRPN encoding small nuclear ribonucleoprotein polypeptide N, rs28915687 in TPO and rs113141381 in YEATS2 encoding YEATS domain-containing 2 may lead to over expression of the respective mRNA transcripts by completely disrupting the MRE so that the miRNA cannot bind (indicated by the higher negative ΔΔG value). Only one SNP, namely rs28562651 in NTRK3 encoding neurotrophic tyrosine kinase receptor type 3, was found to enhance the binding efficiency of hsa-miR-489 to its MRE (indicated by the positive ΔΔG value) thus “over” repressing the translation of NTRK3. Graphical illustrations of each category of MREm-SNPs are provided in Figure 3 (A−C).

Table 1

High-confidence MREm-SNPs in the 3′UTR of autism-implicated genes

Gene symbol	Transcript ID	MREm-SNP ID	Allele change	Putative miRNA affected	miRSVR/total context score	miRNA binding	ΔΔG (kcal/mol)	GERP	phastCons
CSNK1D	NM_139062	rs56124628	G→C	miR-671-5p	−0.692/−0.1	↓	−1.1	3.37	1.0
NTRK3	NM_001007156	rs28562651	T→G	miR-489	−0.422/−0.17	↑	4.1	2.68	1.0
PLCB1	NM_182734	rs28390202	C→T	miR-423-5p	−0.379/−0.4	↓	−2	−0.69	0.991
PRKCB	NM_002738	rs148637931	C→T	miR-194-3p	−0.225/−0.02	↓	−2.2	3.65	0.034
SLC16A3	NM_004207	rs114517919	C→T	miR-550a-5p	−0.208/−0.27	↓	−1.7	1.26	1.0
SNRPN	NM_003097NM_022805NM_022806NM_022807NM_022808	rs112904440	C→G/T	miR-134	−0.734/−0.27	X	−18.8	−3.27	0
TPO	NM_000547NM_001206744	rs1042589	C→G	miR-370	−1.017/−0.44	↓	−3.7	−0.71	0
	NM_001206745NM_175719NM_175721NM_175722	rs28915687	A→T	miR-127-5p	−1.017/−0.15	X	−18.2	−1.38	0
YEATS2	NM_018023	rs113141381	C→T	miR-935	−0.122/−0.29	X	−27.3	1.55	0

Note: MREm-SNP stands for MRE-modulating SNP; decrease and increase in miRNA binding are indicated with ↓ and ↑, respectively, while X indicates complete disruption of miRNA binding. ΔΔG indicates difference in minimum free energy. GERP, genomic evolutionary rate profiling; phastCons, phylogenetic analysis with space/time model for conservation; CSNK1D, casein kinase 1, delta; NTRK3, neurotrophic tyrosine kinase, receptor type 3; PLCB1, phospholipase C beta 1; SLC16A3, solute carrier family 16 member 3; SNRPN, small nuclear ribonucleoprotein polypeptide N; TPO, thyroid peroxidase; YEATS2, YEATS domain-containing 2.

Figure 3

Modulation of existing MREs and creation of new MREs by 3′UTR SNPs in autism-implicated genes The position of SNP within the 3′UTR is indicated by arrow marks. The minimum free energy (MFE) of the miRNA:MRE duplex is enclosed within rectangular boxes. A. The variant allele ‘G’ of rs1042589 decreases the MFE and the binding efficiency of hsa-miR-370 to TPO. B. The variant allele of rs28562651 increases the MFE and binding efficiency of hsa-miR-489 to the 3′UTR of NTRK3. C. Shown is the complementary base pairing between hsa-miR-134 and its MRE in the 3′UTR of SNRPN. The presence of the variant allele ‘G’ completely abrogates the seed paring. D. The variant allele ‘G’ of rs11122391 in the 3′UTR of DISC1 creates a perfect match for the seed region of hsa-miR-137. However, in the presence of the ancestral allele ‘A’, hsa-miR-137 cannot bind the 3′UTR of DISC1. TPO, thyroid peroxidase; NTRK3, neurotrophic tyrosine kinase receptor type 3; SNRPN, small nuclear ribonucleoprotein polypeptide N; DISC1, disrupted in schizophrenia 1.

SNPs in 3′UTR creates new MREs

While concordant cross-prediction of existing MREs and the corresponding binding miRNAs is straight forward, the challenge is to predict those MREs that are created exclusively by the presence of variant alleles at the SNP loci. We followed a logical framework to refine the MREc-SNPs identified with a MAF > 0.1%. First, the SNP involved 3′UTR sequence should not serve as an MRE for any miRNA when the ancestral allele is present, and this was cross verified by (i) the miRanda algorithm embedded within MirSNP, (ii) TargetScan and (iii) RNAhybrid. Second, the introduction of variant alleles at the SNP loci should create a perfect seed match for miRNA(s) and confer favorable hybridization energy for the formation of miRNA:MRE duplex. This was accounted for by analyzing the sequence of each created MRE and its corresponding miRNA for canonical hybridization features in RNAhybrid. The results were manually inspected, and only the MREs that were created for miRNAs expressed in the human brain were considered. A total of 12 SNPs were found to create 35 new MREs across the 3′UTR of 29 transcripts with favorable hybridization energy for 14 different miRNAs (Table 2). These MREc-SNPs may lead to down regulation in the expression of their corresponding genes by providing illegitimate MREs for additional new miRNAs. Of these, rs7701616 created an MRE for hsa-miR-369-3p across 17 transcripts of the PCDHA gene cluster encoding protocadherin alpha. Of note, the SNP rs140716866 in the 3′UTR of PIK3CG encoding phosphotidylinositol-4, 5-bisphosphate 3-kinase, catalytic subunit gamma created a new MRE to which three neuronal miRNAs, hsa-miR-140-3p, hsa-miR-589-3p and hsa-miR-497-3p, can bind, while rs138867267 in the 3′UTR of SCFD2 encoding sec 1 family domain-containing 2 created a new MRE for hsa-miR-33a-5p and hsa-miR-33b-5p. The remaining SNPs created a MRE for one miRNA each and a graphical representation of an example is given in Figure 3D.

Table 2

High-confidence MREc-SNPs in the 3′UTR of autism-implicated genes

Gene symbol	Transcript ID	MREc-SNP ID	Allele change	miRNA with MRE created	ΔΔG (kcal/mol)	GERP	phastCons
DISC1	NM_001164541	rs11122391	A→G	miR-137	18.1	−0.52	0
	NM_001164547	rs61737326	C→T	miR-1255a	16.6	0.31	0
EPHB2	NM_004442, NM_017449	rs2229871	C→T	miR-143-3p	20.3	3.04	0.999
HEPACAM	NM_152722	rs3802904	A→G	miR-1292	29.4	0.86	0
HLA-DRB1	NM_002124	rs9269693	G→C	miR-17-3p	24.9	−0.26	0.017
PCDHA1	NM_018900, NM_031411	rs7701616∗	A→G	miR-369-3p	16.9	3.15	1.0
PIK3CG	NM_002649	rs140716866	A→G	miR-140-3p	18.7	−1.8	0.027
				miR-589-3p	17.1	−1.8	0.027
				miR-497-3p	18.6	−1.8	0.027
RGS7	NM_002924	rs116246787	T→C	miR-143-3p	17.5	2.94	0.757
SCFD2	NM_152540	rs138867267	G→A	miR-33a-5p	20.2	−8.21	0
				miR-33b-5p	18.6	−8.21	0
TNIP2	NM_001161527, NM_024309	rs16843312	G→A	miR-1262	19.2	−1.73	0.001
XIRP1	NM_001198621	rs3732383	G→A	miR-301a-5p	26.2	0.24	0.016
		rs75731397	G→A	miR-185-3p	33.2	−7.54	0

Note: MREc-SNP indicates MRE-creating SNPs; ΔΔG means difference in minimum free energy. ‘∗’ indicates that this SNP also creates MRE for miR-369-3p across other PCDHA family genes, including PCDHA2 (NM_018905), PCDHA3 (NM_018906), PCDHA4 (NM_018907), PCDHA5 (NM_018908), PCDHA6 (NM_018909), PCDHA7 (NM_018910), PCDHA8 (NM_018911), PCDHA9 (NM_014005), PCDHA10 (NM_018901, NM_031860), PCDHA11 (NM_018902), PCDHA12 (NM_01890), PCDHA13 (NM_018904), PCDHAC1 (NM_018898) and PCDHAC2 (NM_018899). GERP, genomic evolutionary rate profiling; phastCons, phylogenetic analysis with space/time model for conservation; DISC1, disrupted in schizophrenia 1; EPHB2, EPH receptor B2; HEPACAM, hepatic and glial cell adhesion molecule; HLA-DRB1, major histocompatibility complex class II, DR beta 1; PCDHA1, protocadherin alpha 1; PIK3CG, phosphoinositide-3-kinase, catalytic subunit gamma; RGS7, regulator of G-protein signaling 7; SCFD2, sec1 family domain-containing 2; TNIP2, TNFAIP3-interacting protein 2; XIRP1, xin actin-binding repeat-containing 1.

Aberrant gene expression and associated mouse phenotypes

As described before, we selected only the genes implicated in autism and miRNAs expressed in the human brain. Therefore, interference in miRNA regulation caused by MRE-SNPs may lead to neuro-pathological effects ranging from low to very high. To better understand the consequences of altered gene expression induced by the MRE-SNPs, we compiled data on the phenotypes produced in mice when the gene homologs were knocked out or targeted. The mouse phenotype data were retrieved from Mouse Genome Informatics (MGI) database [43]. Mice with homozygous null or targeted genes, such as Ntrk3, Plcb1, Prkcb, Disc1, Ephb2 and Rgs7, were shown to exhibit deficits or changes in various neurobiological activities (Table 3). Mouse defective for three of these genes were associated with human disease models. These genes include Plcb1 for Alzheimer’s disease and schizophrenia, Snrpn for Angelman Syndrome and Prader-Willi Syndrome, and Disc1 for schizophrenia and major depressive disorder. On the other hand, the phenotypes of mice defective for Csnk1d, Slc16a3, Tpo, Yeats2, Ephb2, Hepacam, H2-Eb1 (homolog for HLA-DRB1), Pcdha1, Pik3cg, Scfd2, Tnip2 and Xirp1 were not associated with any observable neurobiological deficits (Table S1). However, these genes are associated with the etiology of ASD in humans. Altogether, our results provide in silico evidence for the potential of MRE-SNPs to induce aberrant expression of ASD implicated genes in the human brain that might lead to neurobiological dysfunction and cognitive deficits. The identified SNPs may also add to the genetic heterogeneity besides their role in ASD predisposition and pathogenesis in at least a subpopulation of individuals. These candidate SNPs merit future experimental investigations.

Table 3

ASD candidate genes associated with mouse phenotypes exhibiting neurobiological deficits

Gene symbol	MGI ID	Associated mouse phenotypes	Associated human disease model
Ntrk3	97385	Homozygotes for targeted mutations show a range of phenotypes including postnatal death at 2–21 days, cardiac defects, reduced numbers of dorsal root ganglia neurons and germ cells, abnormal motor coordination and posture and abnormal sensory innervation
Plcb1	97613	Homozygotes for a targeted null mutation exhibit spontaneous seizures and high mortality around 3 weeks of age. Mutant males exhibit sperm with a reduced acrosome reaction rate and fertilizing capacity in vitro and decreased fertility in vivo	Alzheimer’s disease, schizophrenia
Prkcb	97596	Mice homozygous for a null allele exhibit impaired humoral immune responses, altered proliferative responses of B cells to various stimuli, abnormal vascular wound healing and deficits in contextual and cued fear conditioning. ENU-induced mutations lead to impaired T cell-independent IgM responses
Snrpn	98347	Homozygotes for targeted intragenic deletions are phenotypically normal. Deletions that also encompass neighboring genes on the paternal chromosome exhibit growth retardation, hypotonia and high mortality	Angelman syndrome, Prader-Willi syndrome
Disc1	2447658	Homozygotes for a null allele show altered anxiety, synaptic depression, LTP, impulsivity, social investigation, hyperactivity and prepulse inhibition. Homozygotes for a spontaneous allele show altered working memory. Different ENU mutations cause distinct depression and schizophrenia-like profiles	Schizophrenia, major depressive disorder
Ephb2	99611	Mice homozygous for a null allele exhibit abnormal axon guidance, circling, head bobbing and hyperactivity
Rgs7	1346089	Mice homozygous for a hypomorphic allele exhibit reduced exploration in a new environment, impaired glucose tolerance in males and abnormal rod b-wave electrophysiology. Mice homozygous for a knock-out allele exhibit runting, delayed eye opening and transient prolonged b-wave implicit time

Note: Data were compiled from Mouse Genome Informatics (MGI) database. Deficits or changes in neurobiological functions in the mice models are italicized and underlined. Genes associated with human disease models are highlighted in bold. ENU, N-ethyl-N-nitrosourea; Ntrk3, neurotrophic tyrosine kinase receptor type 3; Plcb1, phospholipase C beta 1; Prkcb, protein kinase C beta; Snrpn, small nuclear ribonucleoprotein N; Disc1, disrupted in schizophrenia 1; Ephb2, Eph receptor B2; Hepacam, hepatocyte cell adhesion molecule; Rgs7, regulator of G protein signaling 7.

Discussion

Despite recent advances in genetics and imaging, the neurobiology of ASD remains largely obscure indicating that there is ample space for exploring new hypotheses regarding the etiological basis of the disorder [44]. SNPs are the most common form of genetic variation in the human genome and are known to confer risk and heterogeneity to human disease pathogenesis, clinical course, and response to treatment [45,46]. Owing to its polygenic etiology, SNPs in several genes have been associated with ASD [47]. Apart from the SNPs that occur in protein-coding regions of the gene, regulatory polymorphisms in non-coding regions of the genome might also play a role in driving human phenotypic variation [48,49]. Many autism-implicated genes are under the post-transcriptional regulation of miRNAs, and several human genetic studies have provided some potential connections between miRNA abnormalities and ASD phenotypes. Using multiplex quantitative PCR, Abu-Elneel et al. [19] compared the expression of 466 human miRNAs in postmortem cerebellar cortex tissue samples obtained from 13 ASD individuals and an equal number of non-autistic control individuals. Expression of 28 miRNAs in ASD patients differed significantly from that in the non-autism control set in at least one of the autism samples. In another study [50], microarray profiling identified 9 miRNAs to be differentially expressed in lymphoblastoid cell line samples of ASD patients as compared to the matched controls. Additionally, a number of miRNA expression profiling studies have examined miRNA dysregulation in lymphoblastoid cell lines of ASD patients [51,52]. Recently, our group addressed the biological and functional significance of miRNAs residing at autism-associated copy number variants and their target genes by using an array of computational tools [53]. Given the involvement of miRNAs in ASD pathogenesis, characterization of SNPs that occur in cis sequences, such as MREs, is essential for understanding this peculiar mode of miRNA dysfunction. This mode of miRNA dysfunction does not directly alter miRNA or mRNA abundance, but interferes with their interaction instead. Of note, the first psychiatric disease linked to genetic variation in the miR-189 MRE of SLITRK1 encoding SLIT and NTRK-like family member 1 was the Tourette’s syndrome [25]. Moreover, in a recent study, the promoter and 3′UTR of NLGN4 (neuroligin 4) in a group of 285 ASD cases and 384 controls were sequenced and some 3′UTR variants such as rs3810688, rs5916269 and rs140700235 were found within putative binding sites of miRNAs that might affect NLGN4 expression [54].

Numerous databases [31,55,56] have cataloged several thousands of MRE-SNPs, and identification of functional SNPs amid these MRE-SNPs associated with a complex trait like autism is extremely challenging. We developed a rigorous computational approach based on the latest knowledge of miRNA-mediated post-transcriptional gene regulation. Using this approach, we identified a panel of high-confidence candidate SNPs that may perturb the expression of autism-implicated genes. The MREs predicted to be modulated by SNPs have a good miRSVR score and total context score indicating their functional significance. The MREs created exclusively by the variant alleles of SNPs provide a perfect match for the miRNA seed sequence and favorable hybridization energy, thus the corresponding genes may be “over” repressed from translation. The phastCons and GERP values are also provided to underscore the level of conservation of the SNP loci. Of further interest is the characteristic MAF of these MRE-SNPs (Table S2). For example, some SNPs such as rs28562651 and rs11122391 have a very low MAF in EAs compared to AAs (0.016% in EA vs. 6.77% in AA and 0.15% in EA vs. 24.34% in AA, respectively). The converse is also true as exemplified by rs1042589 (44.37% in EA vs. 23.84% in AA) and rs3732383 (22.4% in EA vs. 4.72% in AA). Moreover, certain SNPs like rs148637931, rs112904440, rs61737326, rs2229871, rs140716866 and rs11624787 were detectable only in AAs. This demonstrates that some SNPs are positively selected in one population and under purifying selection in the other population. The SNPs rs1042589, rs3802904 and rs9269693 are common (MAF > 5%), while a few are uncommon (MAF < 5%) and many others are rare (MAF < 0.5%). Thus a precise estimate of the relative contribution of rare and common variation to ASD risk will require further research [57].

Not only the genes, the miRNAs known to bind the SNP-involved MREs identified in their 3′UTR were also previously known for their association with the pathophysiology of ASD and other neuronal disorders. For example, the microRNAs hsa-miR-671-5p, hsa-miR-489, hsa-miR-589-3p, hsa-miR-497-3p, hsa-miR-33b-5p and hsa-miR-185-3p were reported to reside within CNVs associated consistently with ASD [53]. Sarachana et al. [51] identified hsa-miR-185 and hsa-miR-194 to be up regulated in lymphoblastoid cell lines of autistic individuals. Also, altered levels of hsa-miR-550 and hsa-miR-140 have been implicated in autism [19]. The miRNAs hsa-miR-137, hsa-miR-301a, hsa-miR-489 and hsa-miR-497 have been associated with depression [58]. In addition, all the above mentioned miRNAs except hsa-miR-33b-5p are expressed in the post-mortem brain tissue samples of individuals affected by schizophrenia or bipolar disorder [59]. Overall, our results as well as other studies indicate that miRNAs can have great impact on neuronal function and communication by regulating the expression of their target genes and the presence of MRE-SNPs adds another layer of complexity to the genetic architecture of neuropsychiatric disorders such as ASD, schizophrenia, etc.

Our study has carefully considered the essential features that determine miRNA-mediated gene regulation to better understand the effect of MRE-SNPs and their association with ASD. Despite the requirement of a perfect match between miRNA seed region and MRE [60], it has been shown that G:U wobble in the seed region is acceptable [61,62]. To account for this, we admitted G:U wobble in seed pairing while predicting the miRNA:MRE interactions in RNAhybrid, so that the variant allele of the SNP that simply leads to G:U wobbling without any change in the miRNA:MRE pairing and hybridization energy can be excluded. As the miRNA:mRNA duplex is known to contain self-loops and bulge loops, SNPs in these regions may not have an effect on base pairing or hybridization energy and hence they were omitted. For a subset of SNPs, there was ambiguity in the ancestral allele when cross-checked with NCBI’s dbSNP, and such cases were also excluded. Wherever discrepancies existed in prediction between miRanda, TargetScan and RNAhybrid, the results were discarded. Despite these considerations, a few limitations exist. Human neurodevelopment is a complex and precisely-regulated process, in which nearly 90% of the expressed genes and a set of miRNAs show differential regulation across brain regions and/or time [42,63]. Therefore, spatio-temporal expression of miRNAs and their target genes may confound the effect of these MRE-SNPs. Further, we have focused only on MRE-SNPs, leaving behind the smallest fraction of 3′UTR SNPs that affect 3′UTR splicing. These SNPs are predicted to mediate loss of MREs, mostly through the gain of acceptor splice sites resulting in shortened 3′UTRs [64]. Similar to MRE-SNPs, SNPs in miRNA sequence may also perturb miRNA binding and subsequent regulation of gene expression, or they can affect the processing or maturation of the miRNAs, thus resulting in disease susceptibility [24,65]. Besides RNA–RNA interactions, RNA binding proteins might play a role in miRNA target site recognition as well, either via blocking or facilitating access to the target region [66]. As the list of miRNA-related SNPs and other regulatory variants keeps growing, integrative analysis and experimental validation of their effect on gene expression will be useful for delineating the complex etiology of autism and for deciphering how the different etiological factors act in concert to bring about the various phenotypes. This work points out a unique genetic mechanism that alters the translation of genes in the autistic brain and also provides potential candidates for future experimental validation. Our data suggest that MRE-SNPs may contribute to autism susceptibility, genetic heterogeneity and phenotypic variability. Therefore, MRE-SNPs hold promise for a better understanding of ASD with a consequent improvement in the tools for diagnosis and treatment.

Materials and methods

Autism-implicated gene dataset

The genes implicated in autism were retrieved from the Human Gene Module of the SFARI Gene, which serves as a comprehensive, up-to-date reference for all known human genes associated with ASD [30]. SFARI Gene is a publicly available, curated, web-based, searchable database built by extracting information from the studies on molecular genetics and biology of ASD and can be accessed at http://gene.sfari.org. At the time of this study, 484 genes (Table S3) were listed in the “Human Gene Module” as autism-implicated genes (accessed on April 8, 2013).

Preliminary characterization of MRE-SNPs

The official gene symbols of the autism-implicated genes were uploaded in text format to MirSNP database [31] to retrieve the list of MRE-SNPs. The MirSNP database, available at http://202.38.126.151/hmdd/mirsnp/search/, was developed by integrating mRNA sequence information from NCBI, SNPs from dbSNP135 and miRNAs from miRBase 18.0. MirSNP employs miRanda miRNA target prediction algorithm [34], which incorporates current biological knowledge on target rules, and an up-to-date compendium of mammalian miRNAs. The MRE-SNPs were checked in EVS to obtain information on their MAFs, phastCons scores and GERP scores. The EVS [32] is maintained by NHLBI Exome Sequencing Project and is freely accessible at http://evs.gs.washington.edu/EVS/ with options for downloading the datasets. Only the SNPs with a MAF greater than 0.1% were considered for further analyses.

Identification of functional MREs and MREm-SNPs

The SNP-involved MREs predicted by MirSNP’s miRanda algorithm were subjected to cross-prediction with two additional miRNA target site prediction algorithms namely TargetScanHuman 6.2 [37] and RNAhybrid [38,39]. TargetScan prediction criteria is based on strong pairing in the seed region, thermodynamic stability, number of target sites on the 3′UTR of a given mRNA, sequence context (location of MRE with regard to stop codon and to the poly A tail or presence of AU-rich clusters) and accessibility of the target site to the RISC complex [13,36]. Besides strong/moderate pairing in the seed region and thermodynamic stability, RNAhybrid offers the unique possibility of tuning the search parameters like seed length and the presence of G:U wobble in seed pairing [38,39]. Therefore, this cross-prediction approach covers the multiple features shown to be informative for MRE identification and thus reduces false positives. In addition, the MREs were ranked and selected based on their miRanda-miRSVR score [40] and TargetScan’s total context score [37]. In other words, only the MREs with a “good” miRSVR score of <−0.1 were considered and whenever this score was unavailable, those with a good total context score of ⩽−0.1 were included. Further, the hybridization energy required for the formation of miRNA:MRE duplex was calculated by submitting the sequence of 3′UTR segments containing the MREs and their respective miRNAs in RNAhybrid. Only the duplexes with favorable hybridization energy of ⩾−15 kcal/mol were chosen. The SNP induced mismatch in miRNA:MRE interaction was then analyzed in RNAhybrid by replacing the ancestral allele with the derived allele at the SNP loci of the MRE. The difference in minimum free energy (MFE) of binding before and after introduction of the SNP was computed as ΔΔG (kcal/mol) using the following formula:

In case of SNPs that totally disrupt the formation of miRNA:MRE duplex, the MFE was set to ‘0’ with respect to the derived allele.

Identification of MREc-SNPs

The list of MREc-SNPs with MAF > 0.1% were verified in a step-wise manner as detailed below. Since the MRE is created only after the introduction of the variant allele, the original sequence of the 3′UTR should be devoid of that particular MRE. This was confirmed by the combined use of TargetScan and RNAhybrid. As these MREs do not have an miRSVR score or a total context score, we verified their effect on miRNA binding by checking the sequence of each created MRE and their respective miRNA for canonical hybridization features in RNAhybrid. The ΔΔG was calculated as previously described except that the MFE of the miRNA:MRE duplex was set to ‘0’ with respect to the ancestral allele.

Shortlisting MRE-SNPs based on tissue-specific miRNAs

The neurobiological basis of ASD indicates that the implicated genes mostly function at the synapse and are expressed in the brain. Likewise, the miRNAs that regulate these genes should also be expressed in the brain. Previously, Hu et al. [41] performed both Illumina small RNA sequencing and Agilent miRNA microarray to detect the miRNAs expressed in human brain tissue samples that were obtained from the NICHD Brain and Tissue Bank for Developmental Disorders. Another recent study [42] assessed the spatio-temporal expression of miRNAs in the prefrontal cortex, hippocampus and cerebellum of 18 normal human donor brains spanning different developmental time periods right from infancy through adulthood. In order to achieve consistency and robustness, we utilized the overlapping list of brain expressed human miRNAs from both studies for shortlisting those SNP-involved MREs that serve as their binding sites.

Mapping aberrant expression of candidate genes with mouse phenotypes

To elucidate the phenotypic consequences of aberrant gene expression caused by the candidate MRE-SNPs, we queried their gene symbols in Mouse Genome Informatics (MGI) database [67] using the Genes and Markers Query Form, and retrieved the mouse phenotypes related to targeted knock out or homozygous null mutation of the homologous allele. MGI is a complete catalog of phenotypic mutations in laboratory mouse, and curates data using several different ontologies including the Mammalian Phenotype Ontology [68]. The mouse genotypes that are explicitly used as a model for human diseases are additionally linked with the records of Online Mendelian Inheritance in Man (OMIM). Further, the candidate genes and the miRNAs known to bind the SNP-involved MREs were checked in the scientific literature for any relevance to autism.

Authors’ contributions

AKM and MM conceived and designed the study. VV and MM collected and analyzed the data. MM, VV and AKM wrote the paper. All authors read and approved the final manuscript.

Competing interests

The authors declare no competing financial interests.