Comparative RNA editing in autistic and neurotypical cerebella.

Adenosine-to-inosine (A-to-I) RNA editing is a neurodevelopmentally regulated epigenetic modification shown to modulate complex behavior in animals. Little is known about human A-to-I editing, but it is thought to constitute one of many molecular mechanisms connecting environmental stimuli and behavioral outputs. Thus, comprehensive exploration of A-to-I RNA editing in human brains may shed light on gene-environment interactions underlying complex behavior in health and disease. Synaptic function is a main target of A-to-I editing, which can selectively recode key amino acids in synaptic genes, directly altering synaptic strength and duration in response to environmental signals. Here, we performed a high-resolution survey of synaptic A-to-I RNA editing in a human population, and examined how it varies in autism, a neurodevelopmental disorder in which synaptic abnormalities are a common finding. Using ultra-deep (>1000 × ) sequencing, we quantified the levels of A-to-I editing of 10 synaptic genes in postmortem cerebella from 14 neurotypical and 11 autistic individuals. A high dynamic range of editing levels was detected across individuals and editing sites, from 99.6% to below detection limits. In most sites, the extreme ends of the population editing distributions were individuals with autism. Editing was correlated with isoform usage, clusters of correlated sites were identified, and differential editing patterns examined. Finally, a dysfunctional form of the editing enzyme adenosine deaminase acting on RNA B1 was found more commonly in postmortem cerebella from individuals with autism. These results provide a population-level, high-resolution view of A-to-I RNA editing in human cerebella and suggest that A-to-I editing of synaptic genes may be informative for assessing the epigenetic risk for autism.

Introduction

Site-specific adenosine-to-inosine (A-to-I) RNA base conversions, carried out by adenosine deaminase acting on RNA (ADAR) enzymes, exhibit precise regional specificity in the brain and modulate complex behavior in model organisms[1-3]. A-to-I RNA editing is an efficient means to increase RNA complexity, thereby fine-tuning both gene function and dosage[4-6]. The cellular machinery recognizes inosine as guanosine, so A-to-I editing of codons and splicing signals directly modifies protein-coding gene function[6-11], while editing of microRNAs[12-14] and their binding sites[15] alters gene expression. This is particularly important in the human brain, the single most complex organ in cellular diversity, connectivity, morphogenesis, and responses to environmental stimuli[16].

Synaptic function is a major target of A-to-I editing[17], which can fine-tune neurophysiological properties in response to environmental stimuli[18, 19]. Canonical signaling pathways acting on the editing enzymes link A-to-I RNA editing to environmental cues: ADARB1 function requires inositol hexakisphosphate[20], and the expression of ADAR is interferon-inducible[21]. Several recoding events that directly alter synaptic strength or duration in response to environmental signals have been characterized in rodents[7-11]. For example, mRNAs of the serotonin receptor HTR2C undergo editing in five sites, which dramatically alters its G-protein coupling activity, and hence the relationship between serotonin levels and postsynaptic signal transduction[2, 11]. This editing is regulated by exposure to acute stress and chronic treatment with antidepressants[22]. Another example is the neurodevelopmentally-regulated editing of transcripts encoding the AMPA receptors GRIA2, GRIA3 and GRIA4[23], where arginine to glycine (R/G) recoding of the ligand binding domains leads to faster desensitization recovery[7]. Moreover, glutamine to arginine (Q/R) editing of the transmembrane domains in the kainate receptors GRIK1 and GRIK2 reduces their calcium permeability[8], with varying degrees of editing throughout mouse neurodevelopment[23]. Since 0 to 100% of mRNA molecules can be edited at any given point[24], Q/R and R/G editing of ionotropic glutamate receptors provides an efficient means for fine-tuning the glutamatergic synapse in response to the changing environment.

Little is known about A-to-I RNA editing in humans, but it has been postulated to be one of the molecular mechanisms connecting environmental inputs and behavioral outputs[18, 25]. The increased editing in humans[26] as compared to other animals[27], including nonhuman primates[28], has been proposed to generate molecular complexity that might constitute the basis of higher-order cognition[18]. Therefore, characterizing A-to-I editing in typically and atypically developed individuals may shed light on environment-dependent epigenetic mechanisms central to human neurodevelopment.

Autism Spectrum Disorders (ASD; [MIM 209850]) are highly heritable common neurodevelopmental disorders of complex genetic etiology, characterized by deficits in reciprocal social interaction and repetitive behaviors[29]. Several studies characterize synaptic abnormalities in ASD[29-32]. However, the mechanisms for gene-environment interactions and their contribution to the observed synaptic alterations remain unknown. The number of candidate genes is rapidly increasing[33-35] and a main challenge is to identify the context in which they confer risk to ASD. Recent twin studies suggest that the contribution of environmental factors to ASD is larger than previously thought, with lower monozygotic concordance estimates (77 to 88%), and a surprisingly high dizygotic concordance of 31% (as compared to a sibling recurrence rate of 19%)[36, 37]. Hence, the identification of mechanisms governing gene-environment interactions relevant to ASD could be informative for risk assessment[38]. A-to-I RNA editing is potentially one such mechanism, linking environmental stimuli with synaptic transmission.

Several lines of evidence support an examination of the link between A-to-I editing and ASD. First, model organisms with altered A-to-I editing exhibit maladaptive behaviors characteristic of ASD[2, 39], sometimes with seizures[1, 40] or Prader-Willi-like symptoms[41], both of which are typically detected in 25% of children with ASD[42, 43]. Second, altered editing of the serotonin receptor HTR2C has been detected in a mouse model of autism[44] and in disorders that aggregate in families with ASD[45], including schizophrenia[46] and major depression[47]. Third, a fly model of Fragile × syndrome, the most common single gene cause of ASD[48], was recently shown to exhibit significant editing alterations, via a direct interaction between the Fragile × gene dFMR1 and the editing enzyme dADAR[49]. Finally, independent genomic studies have implicated variants in synaptic genes, the most edited type of genes[17], as a recurring theme in ASD[29-32].

Here we focus on neurodevelopmentally-regulated A-to-I editing that directly alters synaptic function. We precisely quantify and compare the levels of editing across individuals, and characterize the distinct editing landscapes of ten synaptic genes acting in the human cerebellum, contrasting postmortem brains from typically developed individuals with brains isolated from individuals with ASD. We then (i) specifically examine editing patterns in the glutamatergic and serotonergic systems, (ii) correlate editing with isoform usage, and (iii) identify clusters of correlated sites. Importantly, we find that the relative usage of a dysfunctional form of the editing enzyme ADARB1 is significantly higher in ASD.

Materials and Methods

Subjects

Thirty fresh-frozen cerebellar samples of deidentified individuals with nonsyndromic autism and carefully-matched neurotypical individuals were obtained from the National Institute of Child Health and Human Development Brain and Tissue Bank (NICHD BTB) and the Harvard Brain Tissue Resource Center (HBTRC), through the Autism Tissue Program (ATP). To minimize confounding factors, matching was based on age, gender, race, and post-mortem-interval (PMI). Supplementary Table 1 details the 25 samples that passed quality control measures.

Selection of Target Genes

With the accumulation of multidisciplinary studies highlighting synaptic alterations in ASD[29, 31, 32, 34, 50–54], we chose to focus on neurodevelopmentally-regulated A-to-I recoding that alters synaptic function. All synaptic genes shown to undergo developmentally dependent A-to-I recoding by Wahlstedt et al.[23] were included in this study (Supplementary Table 2).

Molecular Methods

RNA isolation and quality assurance

Tissue samples were disrupted and RNA isolated in triplicate (mirVana™ miRNA Isolation Kit, Life Technologies, Grand Island, NY, USA), followed by DNase I treatment (DNA-free™, Life Technologies). Samples that seemed intact on 1% agarose gels were quantified on the ND-1000 Spectrophotometer (NanoDrop, Wilmington, DE, USA) and their integrity analyzed using the Agilent 2100 Bioanalyzer Eukaryote Total RNA Nano Series II (Agilent Technologies, Santa Clara, CA, USA). Only samples with RNA Integrity Number (RIN) >7 were included in the study, minimizing post mortem degradation artifacts[55]. The mean RIN of the samples used was 8.0 (Supplementary Table 1).

454 cDNA library preparation

A stringent PCR-based library preparation protocol was designed to ensure unbiased templates for multiplexed 454 sequencing, as detailed in Supplementary Information and Supplementary Figure 1.

454 sequencing

Bidirectional GS FLX sequencing was performed as described[56] by the 454 Life Sciences Sequencing Center (Branford, CT, USA). ASD and neurotypical pools were sequenced on opposite sides of a 2-region PicoTiterPlate. 457,104 reads were obtained, containing 85,709,299 high quality bases (Supplementary Figure 2).

DNA isolation

About 2g of the same cerebellar tissue samples were disrupted with a mortar and pestle on dry ice, and genomic DNA (gDNA) was extracted using the QIAamp DNA Mini Kit (QIAGEN, Valencia, CA, USA).

Genomic DNA genotyping

Sequenom iPLEX genotyping (Sequenom, Inc, San Diego, CA, USA) was done at Boston Children’s Hospital's SNP Genotyping Facility. Twenty-five SNPs were genotyped in four populations. Every sample was genotyped in triplicate, and all triplicates had to agree on the genotype to be considered successful. All passing SNPs had a genotype success rate of >96.7%.

Method validation

cDNA and gDNA from all samples were independently analyzed by parallel capture and Illumina sequencing, as detailed in Supplementary Information. On average, 3420 cDNA reads and 627 gDNA reads were obtained for every sample at each site tested. The results of the two editing detection methods were tightly correlated, with an average Pearson’s r of 0.923 (mean p= 0.002) across six sites. See Supplementary Table 4 for details.

Relative ADARB1 isoform usage analysis

Semiquantitative RT-PCR and Sanger sequencing were used to measure the relative frequencies of the normal and dysfunctional (NR_027672) ADARB1 isoforms, as detailed in Supplementary Information. Custom TaqMan® Gene Expression Assays were used to quantitate the abundance of the dysfunctional ADARB1 isoform relative to GAPDHas detailed in Supplementary Information.

Correlation between editing and splicing of FLNA

When we amplified across the FLNA Q/R site using one primer in exon 43 and one in exon 44, an unexpected 720bp band of varying brightness appeared in all samples. Bidirectional Sanger sequencing revealed that it is a result of intron 43 retention. Its relative abundance was quantitated using QIAGEN’s QIAxcel electrophoresis system, and correlated with Q/R editing levels detected by 454 sequencing.

Gene expression profiling

Affymetrix Exon 1.0 ST arrays (Santa Clara, CA, USA) were used to measure global gene expression of each RNA sample. Following quantile normalization, gene expression levels were calculated using the Probe Log Iterative ER (PLIER) algorithm and differences in Gene Core expression were determined using the Mann-Whitney U test.

RNA-seq

Ten of the surveyed samples were subject to high-throughput RNA sequencing on Applied Biosystem’s SOLiD platform (Life Technologies). Poly(A)+ RNA was captured (Oligotex mRNA Mini Kit, QIAGEN), heat-fragmented (95°C for 20 minutes), and prepared for sequencing (Small RNA Expression Kit, Life Technologies) that yielded 114M reads. Life Technologies’ WT Analysis software package was used to align the reads to the reference genome and expressed sequence tags, identifying 33.4M uniquely mapped reads.

Computational Methods

454 variant calling and quantification

The 454 GS Amplicon Variant Analyzer software was used to deconvolute samples, align reads to reference sequences, and call variants and haplotypes based on bidirectional sequence changes from the reference.

Beta-binomial modeling of editing levels

For each sample in each site, the posterior editing density, f(editing), was a beta distribution with parameters α=1+number of reads with G mapped to the site, and β=1+number of reads with A mapped to the site (Supplementary Figure 3).

Differential editing analyses

Kolmogorov-Smirnov (KS) tests were used to assess the significance of differences in continuous measurements, such as those of editing levels between neurotypical individuals and individuals with ASD. To compare discrete distributions, such as the count of GRIK2 or HTR2C isoforms resulting from combinatorial RNA editing, Pearson’s chi square test was used. The total number of reads belonging to each of the possible editing isoforms were summed across each group and only isoforms supported by > 5 reads in both groups were included in the analysis (following this test’s assumptions).

Correlation between dysfunctional ADARB1 Isoform frequency and overall editing levels

For each sample i, RF is the relative frequency of the dysfunctional ADARB1 isoform. A sample’s overall editing level, E, is the sum of its standardized editing levels across all sites, ∑Z, where Z is the standardized editing level of sample i at site j. Pearson’s correlation was calculated between these two metrics, RF and E, across all samples, to measure their linear dependence.

Association between editing and splicing of AMPA Receptors

Two-by-two contingency tables were used to summarize the relationships between isoform selection and editing (Supplementary Table 5), and Fisher’s exact test was used to determine the significance of the association between them.

Correlations among editing sites

The Pearson correlation coefficient was calculated to quantify the linear relationships between editing at different sites, among all individuals. Biclustering was then used to identify modules of tightly correlated sites, with the EXPANDER software[57].

Independence of editing at neighboring sites

To summarize the relationships between editing at GRIA2 Q/R and Q/R+4 among all individuals, a two-by-two contingency table is shown in Supplementary Table 5. Fisher’s exact test was used for power and significance calculations.

Multiple testing correction

All p-values were Benjamini-Hochberg corrected for multiple testing[58] to ensure that the false discovery rate of this entire study is below 0.05.

Results

Precise Multiplex Quantitation of A-to-I RNA Editing in Human Cerebella

For a high-resolution view of A-to-I RNA editing in a human brain population, ultra-deep (>1000×) 454 cDNA sequencing and genomic DNA (gDNA) genotyping were used to detect A/G mixtures on cDNA mapping to homozygous genomic A/A (Supplementary Figure 1). This study focused on all ten synaptic genes shown in mouse to undergo neurodevelopmentally-regulated A-to-I recoding[23] that results in well-characterized neurophysiological alterations[7–9, 59–61] (Supplementary Table 2). Editing was measured in postmortem cerebellum, one of several brain regions implicated in ASD[62] by both imaging[63] and autopsy[64] studies. Cerebellar tissue samples were obtained from neurotypical individuals and individuals with non-syndromic ASD, matched by gender, age, race, and postmortem interval (Supplementary Table 1). Pooled, barcoded, bidirectional 454 cDNA sequencing yielded on average 1344 reads of 212bp per individual per amplicon (Supplementary Figure 2). gDNA was genotyped at 18 well-characterized editing sites with known functional consequences and 7 positions aligned to a mixture of guanosines and adenosines on cDNA. A-to-I edited sites were identified by the presence of a bidirectional cDNA A/G mixture at homozygous A/A gDNA positions. Editing levels were modeled by a beta-binomial distribution, resulting in a posterior editing density for each individual at each site. This model considers both the fraction of edited reads and the total number of reads covering a site, to produce a distribution describing the level of editing and our confidence in that measurement (Supplementary Figure 3).

This approach provided a high confidence dataset of synaptic A-to-I RNA editing from individuals with ASD and matched neurotypical individuals, with an average 95% confidence interval length of 0.038. For independent validation, the same cDNA and gDNA were analyzed by parallel capture and >3000× Illumina sequencing in six sites among all individuals (Supplementary Information). A tight correlation between the two editing detection methods is shown in Supplementary Figure 4 and Supplementary Table 4, with an average correlation coefficient of 0.923.

High Dynamic Range of Editing across Sites and Individuals

The editing levels of 25 sites were robustly measured and found to range from 2.5% to 99.6% (Figure 1A and Supplementary Figure 5). Unexpectedly broad distributions of editing levels were found across neurotypical individuals and carefully-matched individuals with ASD in those sites previously shown to be neurodevelopmentally regulated in mice[23] (average coefficient of variation 54.8%). Only one site, GRIA2 Q/R, showed similar levels of editing among all individuals (CV=0.5%), consistent with reports that editing of this site is essential and unchanged throughout development[23, 65]. No relationships between editing levels and age, gender, race, postmortem interval, or RNA integrity were detected (Supplementary Figure 6).

Figure 1

Individuals with ASD at the extremes of the population editing distributions in 20 of 25 sites. (A) High dynamic range of mean editing across sites, from 99.6% to below detection limit. Standard deviation bars, shown in purple, indicate the large variability in neurodevelopmentally-regulated editing among carefully matched individuals. (B) The individual standardized deviations from the mean editing level across all sites show that the extreme of the population editing distributions tend to be individuals with ASD (orange triangles), and they are the major contributor to the large variability shown in (A). Horizontal lines denote the standardized 95% confidence interval of an individual’s posterior editing distribution. Synaptic editing levels at least two standard deviations away from the mean are highly informative of ASD, with a positive predictive value of 78%.

By inspection of the individual posterior editing densities across all sites (Supplementary Figure 5), we noticed that the extremes of most editing distributions are individuals with ASD. To quantify this observation, the point estimates of the editing levels of each individual at each site were transformed to Z scores and those ≥2 or ≤−2 were considered extreme, representing editing levels that are at least two standard deviations away from the mean. In 20 of 25 sites, individuals with ASD were at the extreme of the spectrum of editing seen for that site (Figure 1B). Having extreme synaptic editing levels is highly informative of ASD, with a positive predictive value p(ASD | |Zediting | > 2) of 0.78. Outliers at different sites are different individuals with ASD (Supplementary Figures 5 and 7), and editing of more than three standard deviations away from the mean was specific to individuals with ASD (Supplementary Figure 7). The overall editing variance in ASD was more than two-fold that of neurotypical individuals (ASD variance=0.58, median neurotypical variance of equally sized subsamples=0.24, Brown-Forsythe p=5.6e-3). The major contributors to this increased variance are 30% of the individuals with ASD (Supplementary Figure 8).

Patterns of Editing in the Glutamatergic and Serotonergic Systems

Next, editing patterns were specifically examined in the glutamatergic and serotonergic systems, which have been shown to be implicated in ASD[29, 66]. GRIK2, whose genomic locus has been repeatedly linked and associated with ASD[66], is edited in three sites, leading to the formation of eight protein isoforms that differ in their calcium permeability[23]. Differential GRIK2 editing between individuals with ASD and neurotypical individuals is depicted in Figure 2, both on a single site and isoform-wide levels. Kolmogorov-Smirnov (KS) tests were used to assess the significance of differences in continuous measurements, such as those of editing levels between neurotypical individuals and individuals with ASD. To compare discrete distributions, such as the count of GRIK2 or HTR2C isoforms resulting from combinatorial RNA editing, Pearson’s chi square test was used. The relative frequencies of GRIK2 isoforms were significantly different between individuals with ASD and neurotypical individuals (p<1e-4, chi-square test). Differences were also detected in the editing of GRIA4 R/G (p<3.6e-2, KS testSupplementary Figure 9). Other comparisons of glutamate receptors did not reach statistical significance, as this study is not powered to detect small effect sizes (Supplementary Figure 10).

Figure 2

Differences in editing-mediated GRIK2 isoforms between neurotypical individuals (blue) and individuals with ASD (orange). GRIK2 is edited at three sites: I/V, Y/C (both part of TM1) and Q/R (TM2), altering the receptor’s calcium permeability. Editing at these sites was found to be tightly correlated (Figure 5B). (A) Large variability in GRIK2 editing at each site among all individuals. The extremes of the observed editing spectra are individuals with ASD. Consult Supplementary Figure 5 for more details. (B) Combinatorial editing of the three sites results in eight isoforms, with distinct molecular properties. Significant group differences in GRIK2 isoform distributions between individuals with ASD and neurotypical individuals were identified (p<1e-4, chi-square test). (C) The strongest differences are in IYQ and VCQ isoforms, which are ~1.5-fold over-represented in ASD. Shown are differences in the cumulative density functions (CDFs) of the VCR, ICR, IYQ, and VCQ isoforms.

Editing was also examined in the serotonin receptor HTR2C, which is targeted by the only FDA-approved drugs to treat autistic symptoms, risperidone and aripiprazole. This receptor undergoes RNA editing at five sites, which dramatically alters its G-protein coupling activity. HTR2C editing may create up to 24 different protein isoforms, characterized by distinct molecular and behavioral phenotypes[2, 6]. Comparing the relative isoform frequencies of individuals with ASD and neurotypical individuals revealed significant differences (p<1e-4, chi square test), the strongest in IDV, MNV and ISI isoforms - all under-represented in ASD (Figure 3). However, the variance within ASD is much larger than that between neurotypical and ASD individuals, suggesting that this finding should be considered with caution and revisited with a larger sample size.

Figure 3

Differences in the relative abundance of HTR2C isoforms between individuals with ASD and neurotypical individuals (p<1e-4, chi-square test). 21 of the possible 24 isoforms created by combinatorial editing of five HTR2C sites were detected. The greatest differences were in the IDV, MNV and ISI isoforms, all under-represented in ASD by >2-fold.

Differential Relative ADARB1 Isoform Usage in Individuals with ASD

The sites examined here are predominantly edited by ADARB1 (Supplementary Table 2). To explore the regulatory basis of the observed editing differences, ADARB1 expression and relative isoform usage were examined. In rodents, the function of ADARB1 has been shown to be developmentally regulated via alternative splicing that creates an inactive protein, while its overall expression remains fairly constant throughout development[23, 67]. In humans, a dysfunctional ADARB1 isoform (ADAR2g, NR_027672), resulting from alternative skipping of the exon harboring two double stranded RNA binding domains (dsRBDs), has been found to constitute about 20% of adult cerebellar ADARB1 mRNA[68]. As such, the presence of this dsRBDs-encoding exon was assayed by semi-quantitative RT-PCR in all samples. While individuals with ASD and neurotypical individuals showed similar levels of overall ADARB1 expression (Supplementary Figure 11), the relative usage of the inactive ADARB1 isoform was significantly more common in individuals with ASD (Figure 4 and Supplementary Figure 12, p=3.0e-3, Mann-Whitney U test). Exon skipping was confirmed by bidirectional sequencing in all samples. Ten of the surveyed samples were also subject to RNA-seq, and the relative ADARB1 exon usage determined by RNA-seq correlated with the semi-quantitative RT-PCR results (Supplementary Information). Additionally, an increased relative usage of the dysfunctional ADARB1 isoform in ASD was observed in an independent cohort of nine cerebellar tissue samples (Supplementary Information). The relative frequency of the dysfunctional isoform is correlated with overall editing levels (Pearson’s r= 0.48, p=3.2e-2), as measured by the sum of the standardized editing scores across all sites. Finally, the abundance of the inactive ADARB1 isoform was quantitated in all samples using TaqMan® Gene Expression Assays (Supplementary Information). The dysfunctional form was over-expressed in individuals with ASD (p=7.3e-3, Mann-Whitney U test), with a significant correlation between semi-quantitative and quantitative PCR results (Pearson’s p= 2.0e-5).

Figure 4

Differences in the relative usage of ADARB1 inactive form between individuals with ASD and neurotypical individuals (p=3.0e-3, Mann-Whitney U test). This isoform lacks the enzyme’s dsRBDs and has been shown to be untranslated[68]. ADARB1 is the predominant editing enzyme of the sites analyzed in this study. (A) Semiquantitative RT-PCR was used to measure the relative frequencies of functional (1057 bp product) and dysfunctional (NR_027672, 122 bp product) ADARB1 isoforms, using primers (black triangles) which detect both. (B) The mean relative frequency of the dysfunctional ADARB1 isoform was 46.3% ± 5.6% in ASD and 26% ± 4.5% in neurotypical individuals, similar to that previously reported in neurotypical human cerebellum[68]. All 25 samples included in the RNA editing survey were also part of this analysis. The relative frequency of the dysfunctional isoform is correlated with overall editing levels (Pearson’s p=3.2e-2). Additionally, the abundance of the dysfunctional form as quantitated by TaqMan® Gene Expression Assays was higher in individuals with ASD (p=7.3e-3, Mann-Whitney U test, Supplementary Figure 12), with a significant correlation between quantitative and semi-quantitative PCR results (Pearson’s p= 2.0e-5).

Tight Relationships between Recoding and Splicing in Three AMPA Receptors and Filamin-A

The long 454 reads (Supplementary Figure 2) were used to study regulatory characteristics of human RNA editing, including the previously hypothesized[6; 69] relationships between recoding and alternative splicing selection in three AMPA receptors (GRIA2, GRIA3, and GRIA4) and filamin-A (FLNA). GRIA2, GRIA3, and GRIA4 each contain two mutually exclusive exons that modulate desensitization kinetics, termed flip and flop[. Editing at the 3’ end of the exon immediately preceding the flip/flop module was strongly associated with the flop isoform in GRIA4 (OR=101.9 [95% CI 87.8–118.2]), GRIA3 (OR= 30.6 [26.8–34.8]), and to a lesser extent GRIA2 (OR= 2.6 [2.4–2.8]) (All p-values <1e-300, Fisher’s exact test, Figure 5A and Supplementary Table 5). Consistent with the detection of differential editing at GRIA4 R/G between individuals with ASD and neurotypical individuals, differential GRIA4 isoform usage was also identified (p<3.6e-2, KS test, Supplementary Figure 13).

Figure 5

Relationships between editing and splicing, and among sites. (A) Editing at the 3’ end of the exon immediately preceding the flip/flop module is strongly associated with the flop isoform in GRIA4 (OR=101.9 [95% CI 87.8–118.2]), GRIA3 (OR= 30.6 [26.8–34.8]), and to a lesser extent GRIA2 (OR= 2.6 [2.4–2.8]) (All p-values <1e-300, Fisher’s exact test). (B) The linear correlation of editing levels across sites was measured among all individuals and biclustered to reveal modules of tightly correlated sites, suggesting coregulation by the same factors. The tightest module contains 6 sites in 4 genes, demonstrating highly correlated editing across these sites among all individuals. Consult Supplementary Figure 5 for broader perspective. (C) Three neighboring sites at a 4 bp distance from one another undergo independent editing (p=1, power>0.97, Fisher’s exact test), with GRIA2 Q/R −4 fully unedited, GRIA2 Q/R fully edited, and GRIA2 Q/R+4 variably edited across individuals, between 4–17%.

Another relationship between recoding and splicing was examined in FLNA, where editing at the 3’ end of exon 43 was found to be significantly correlated with intron 43 retention (Pearson’s p<6e-3, Supplementary Figure 14). This intron retention introduces a frameshift and is not included in any human Refseq FLNA isoform. Alignment of the retained intron sequence to expressed sequence tags and short transcriptome reads revealed its expression in several tissues from multiple organisms. Thus it likely represents a conserved splicing event that possibly regulates FLNA dosage. All retained introns were found to contain an unknown A-to-I editing site at their center, verified by gDNA genotyping. Whether this editing is a cause or effect of the intron retention remains to be elucidated.

cis and trans Relationships across Editing Sites

To learn more about the regulation of human RNA editing, clusters of tightly correlated editing sites across all individuals were identified. First, Pearson’s correlation was calculated to quantify the linear dependencies between editing at different sites, among all individuals. Biclustering was then used to identify modules of tightly correlated sites with an average correlation coefficient > 0.7 between all pairs of sites in the cluster, across all samples (Supplementary Table 6). The cluster with the strongest correlation contained 6 sites in 4 genes (GRIA2 R/G, CYFIP2 K/E, GRIA3 R/G-1, and GRIK2 I/V, Q/R and Y/C) (Pearson’s r=0.8, p<1e-4), suggesting that these sites may be co-regulated (Figure 5B). Some neighboring sites, including GRIA2 Q/R and Q/R+4, showed independent editing (Fisher’s Exact p=1, Power>0.97, Figure 5C and Supplementary Table 5), suggesting that their editing is either spatiotemporally distinct and/or carried out by different complexes.

Discussion

This study represents an initial examination in ASD of A-to-I RNA editing, a form of gene-environment interaction that fine-tunes synaptic function in response to environmental stimuli. Human genes that undergo editing are more often involved in processes that are affected in neurodevelopmental disorders[17, 18], and A-to-I editing is overall increased in the human lineage[26-28]. While little is known about the impact of RNA editing in humans, there is compelling evidence supporting a key role for A-to-I editing in modulating complex behavior in animals[1–3, 41, 44, 70]. Mice and flies with altered editing levels recapitulate several behavioral homologues to human ASD[1, 41, 44]. Gene-specific studies in animals demonstrate the causal relationship between editing levels and specific maladaptive behaviors[1, 10], and between specific environmental exposures and editing levels[71]. Therefore, the results presented here shed light on the role of an epigenetic mechanism that connects environmental signals and downstream behavioral outputs, integrating genetic and environmental information.

DNA sequence variation in a number of different loci is strongly associated with ASD, but no individual locus is altered in more than 2% of cases[66]. In contrast, 30% of the individuals with ASD examined here showed extreme levels of synaptic RNA editing, suggesting that A-to-I editing may be both a marker for and mechanism of ASD. It has been recently shown that typical synaptic protein synthesis occurs within an optimal range and deviations in either direction can lead to cognitive impairments[72]. This study suggests that a similar inverse-U relationship might exist between synaptic A-to-I editing levels and the ASD phenotype. In ASD, it now appears that increased variance can be found in different levels of function, ranging from synaptic A-I editing and protein synthesis[72], through mitochondrial[73] and immune[74] function, to cellular and systems neuroanatomy[75]. A plausible generalization of these findings is that many biological processes exhibit increased variance in ASD as a result of many etiologies disrupting ASD function and homeostasis. We speculate that altered A-to-I editing could act as a common compensatory mechanism for the wide range of synaptic abnormalities in ASD, as affected neurons try to maintain synaptic homeostasis. Alternatively, altered editing could be a direct consequence of ASD-causing mutations. In either case, larger sets of genes, brain regions, and individuals should be examined to understand the potential contribution of A-to-I editing to ASD.