Editing of the pre-mRNA for the serotonin receptor 2C (5-HT2CR) by site-specific adenosine deamination (A-to-I pre-mRNA editing) substantially increases the functional plasticity of this key neurotransmitter receptor and is thought to contribute to homeostatic mechanisms in neurons. 5-HT2CR mRNA editing generates up to 24 different receptor isoforms. The extent of editing correlates with 5-HT2CR functional activity: more highly edited isoforms exhibit the least function. Altered 5-HT2CR editing has been reported in postmortem brains of suicide victims. We report a comparative analysis of the connections among 5-HT2CR editing, genome-wide gene expression and DNA methylation in suicide victims, individuals with major depressive disorder and non-psychiatric controls. The results confirm previous findings of an overrepresentation of highly edited mRNA variants (which encode hypoactive 5-HT2CR receptors) in the brains of suicide victims. A large set of genes for which the expression level is associated with editing was detected. This signature set of editing-associated genes is significantly enriched for genes that are involved in synaptic transmission, genes that are preferentially expressed in neurons, and genes whose expression is correlated with the level of DNA methylation. Notably, we report that the link between 5-HT2CR editing and gene expression is disrupted in suicide victims. The results suggest that the postulated homeostatic function of 5-HT2CR editing is dysregulated in individuals who committed suicide.
The serotonin receptor 2C (5-HT2CR) is abundantly expressed in neurons of the central nervous system (CNS), where it modulates diverse processes that affect sleep, sexual activity, feeding behavior, anxiety and depressive state (reviewed in 1,2). The 5-HT2CR belongs to the G-protein-coupled receptor (GPCR) superfamily and is involved in multiple signal transduction pathways through coupling to different G proteins, which allows 5-HT2CR to mediate intracellular responses via various effector molecules (reviewed in 3,4). Several of these second messengers have been reported to influence the activity of ion channels, and thus affect the excitability of the 5-HT2CR-expressing neurons (5).
The 5-HT2CR shows remarkable functional variability that derives from editing of its pre-mRNA by site-specific adenosine deamination (A-to-I pre-mRNA editing) (4), which is catalyzed by RNA-specific adenosine deaminases (ADAR1 and ADAR2) (reviewed in 6,7). The A-to-I editing most frequently occurs in non-coding Alu repetitive elements (8–10). Although the biological relevance of this editing remains uncertain, its overall level has been reported to be dramatically higher in primates than in other mammals, and moreover, in humans compared with non-human primates, thus suggesting a possible contribution of editing to the development of higher-order cognitive processes (8,11,12). Site-specific mRNA editing resulting in re-coding has been reliably confirmed in only a few transcripts, most of which encode proteins that are expressed in the CNS (reviewed in 13,14). Notably, this type of editing occurs at functionally critical positions in proteins involved in synaptic transmission, including 5-HT2CR, many non-NMDA (AMPA and kainate) glutamate receptor ion channel subunits, α3 subunit of GABAA receptor ion channel, and voltage-gated potassium (Kv1.1.) and calcium (Cav1.3) channels. The protein isoforms resulting from mRNA editing profoundly influence neuronal excitability and signal transduction. Hence, it has been hypothesized that RNA editing provides efficient means for fine-tuning neurophysiological properties in response to environmental stimuli (15). In line with this hypothesis, recent work in cultured neurons and hippocampal slices has demonstrated that pharmacologically induced changes in neuronal activity can impact the extent of the ADAR-dependent editing of specific target mRNAs (16,17).
The 5-HT2CR can be edited at up to five closely-spaced (within a 15 nucleotide sequence) adenine positions that have been denoted A, B, E (also known as C′), C, and D sites. Because inosine is read as guanosine during translation (6), editing can alter codons for three amino acids in the second intracellular loop of the receptor (18,19), a region involved in the coupling with G-proteins (20). Combinatorial editing at the five positions can generate up to 32 mRNA variants encoding 24 different receptor isoforms (sites A and B as well as sites E and C are in the same codons). The extent of editing is inversely correlated with the 5-HT2CR function such that the more highly edited isoforms are less active than less extensively edited ones (reviewed in 4). The unedited Ile156-Asn158-Ile160 (INI) isoform possesses considerable constitutive and agonist-stimulated activity. In contrast, when the mRNA is edited, the coupling of 5-HT2CR to G-proteins and its affinity for serotonin are drastically reduced. Specifically, experiments in heterologous expression systems have shown that, compared with the INI isoform, the VSV (Val156-Ser158-Val160) isoform, which is the most common isoform in the human prefrontal cortex (21–23) and can be encoded by two different mRNA variants (ACD or ABCD), shows 4–5-fold decreased coupling to signaling pathways and 4-fold reduced constitutive activity (18,24).
Studies in animals have revealed variations in 5-HT2CR editing patterns among different strains of laboratory mice, indicating that editing can be genetically regulated (25). Experiments in different animal models have also suggested that 5-HT2CR editing is influenced by exposure to acute stress and chronic treatment with antidepressants (26) as well as by other pharmacological treatments (reviewed in 4, although see 27). A recent study of spinal cord injury (SCI) in a rat model has established a link between an editing-mediated increase in the spinal cord 5-HT2CR activity and the recovery of the spinal motoneuron excitability after SCI (28). Taken together, the diverse nature of 5-HT2CR signaling, its involvement in the regulation of complex behaviors, and its remarkable capacity to vary/adjust its activity via editing based on genetic status or in response to environmental impact suggest that 5-HT2CR editing is a homeostatic mechanism that allows 5-HT2CR-expressing neurons to stably maintain their functional characteristics.
Postmortem studies in several laboratories have reported an association between 5-HT2CR editing and suicide (22,23,29–31). Specifically, our previous findings indicate that in the three major psychiatric diseases, major depressive disorder (MDD), bipolar disorder (BPD) and schizophrenia (SZ), suicide is consistently associated with higher level of 5-HT2CR editing (and accordingly lower receptor activity) in the prefrontal cortex (PFC) (22,23,29–31). Suicide is linked to numerous, diverse risk factors. Among the strongest associations is the presence of psychopathology, most notably MDD, BPD, SZ or substance use disorder (32). Genetic epidemiological data provide strong evidence for a heritable component of suicidal behavior which is in part independent from the familial clustering of mental disorders (33). Nevertheless, environmental variables explain a substantial fraction of the total variance in suicidal behavior (34). Recent studies have increasingly emphasized the role of interactions between genetic vulnerability and environmental impact in suicide, suggesting that environmental stressors could influence subsequent behavior by altering gene function through epigenetic re-programming (35,36). Epigenetic processes link environmental factors to gene function via chromatin modifications (i.e. DNA methylation and histone modifications) that regulate gene expression. Recent reports describe suicide-associated alterations of DNA methylation in the promoter regions of several genes and the effects of these alterations on gene expression (reviewed in 37). These findings imply that epigenetic modifications form an interface through which suicide-associated environmental impacts affect gene expression. Little is known, however, about the biological mechanisms which translate gene-environmental interactions into modulation of neuronal functions that ultimately govern complex behaviors, e.g. suicide. The A-to-I editing in general, and 5-HT2CR editing in particular, is a potentially important mechanism of this type as it could be involved in homeostatic regulation of neurons by linking environmental stimuli to signal transduction and neuronal excitability.
Here we investigated (i) whether 5-HT2CR editing is linked to pathways that are important for neuronal activity and (ii) whether this link is altered in the brain of people that committed suicide. We first quantified and compared the levels of 5-HT2CR editing in the PFC between individuals who died of suicide or non-suicide causes using autopsy brain specimens obtained from two independent cohorts. The results confirm the overrepresentation of highly edited mRNA variants in suicide brains that was previously observed in different studies and populations. We then investigated the relationship between 5-HT2CR editing and gene expression in both study cohorts and detected numerous genes whose expression level was associated with editing. The set of these editing-associated signature genes was significantly enriched for genes that are involved in synaptic transmission, genes that are preferentially expressed in neurons, and genes whose expression is associated with the level of DNA methylation. Finally, we found that the link between 5-HT2CR editing and gene expression was (nearly) eliminated in suicide victims.
Overall analysis of RNA editing
We used massive parallel sequencing (MPS) to investigate 5-HT2CR editing in the brain samples obtained from two cohorts described in Materials and Methods, Supplementary Methods and Supplementary Material, Tables S1 and S2. Cohort 1 consisted of subjects who died of suicide (N = 22) or non-suicide-related causes (control group, NC; N = 29). Cohort 2 consisted of patients with a history of MDD who died of suicide (MDDSU, N = 10) or of non-suicide causes (MDDNon-SU; N = 24) and control individuals without any psychiatric diagnoses (DX) who died of non-suicide causes (NC; N = 31). The samples were dissected from sub-regions of the PFC–orbital frontal cortex (OFC; Brodmann area 11) for Cohort 1 and the subgenual anterior cingulate cortex (ACC; Brodmann area 25) for Cohort 2.
On average, 331 924 ± 3494 and 452 450 ± 44 881 (mean ± SE) MPS reads were obtained for each subject in Cohorts 1 and 2, respectively. The high sequencing depth ensured highly reliable detection of all of the 32 5-HT2CR mRNA editing variants, including those with abundance of <0.001% (38).
We detected only small (<1.3%), although significant for some variants, differences in the variant frequencies between the cohorts (Fig. 1, Supplementary Material, Table S3 and Fig. S1). Accordingly, the patterns of 5-HT2CR variant expression were highly similar between the cohorts and were also consistent with the data obtained in our previous studies that had been performed in yet another PFC sub-region, dorsolateral PFC (DLPFC) (22). In particular, only 3 of the 32 mRNA variants (editing combinations ABCD, ABD and NONE) were observed at ∼10% or greater frequency; ABCD was the most highly expressed variant in all three areas of the PFC.
In the principal component analysis (PCA) of all 5-HT2CR editing variants, the ABCD (edited at sites A, B, C, D) and NONE (non-edited variant) aligned along the first principal component which by itself explained 62.7 and 62.6% of the total data variance in Cohorts 1 and 2, respectively (Supplementary Material, Fig. S2). Therefore, in all further analyses, the level of editing in each individual was represented by the frequency of the ABCD variant (hereinafter ‘ABCD’).
ABCD is increased in the brains of suicide patients
In our previous studies in the DLPFC of BPD and SZ patients, we detected more ABCD in suicides than in non-suicides with the same diagnosis but detected no differences between the non-suicides with BPD or SZ and psychiatrically normal controls (22). In our subsequent study of MDD patients, we found significantly less editing in MDDNon-SU than in MDDSU, whereas the MDDSU did not significantly differ from psychiatrically normal controls (39). Thus, both studies have shown that, compared with non-suicides with the same diagnosis, suicides express higher levels of the ABCD variant that encodes a hypoactive 5-HT2CR protein isoform, suggesting that 5-HT2CR editing (and the receptor function) is altered in the PFC of suicide victims regardless of their underlying psychiatric illness.
Here we first asked if our previous findings could be validated in different cohorts and/or different areas of the PFC. In Cohort 1, we detected a higher ABCD frequency in the suicide compared with the control group (P = 0.0217 by Wilcoxon test) (Fig. 2A). In Cohort 2, we detected a higher ABCD frequency in the MDDSU patients compared with MDDNon-SU individuals [P = 0.0216 by Tukey's honest significant difference (HSD) post hoc test]; there was also more ABCD in the MDDSU compared with the NC group but this comparison was not statistically significant (P = 0.3326 by Tukey's HSD post hoc test).
Using logistic regression analysis, we then compared ABCD frequencies between suicides and non-suicides, adjusting for covariates [i.e. age, sex, pH, postmortem interval (PMI) as well as DX in Cohort 2]. We detected higher ABCD in suicides for both cohorts (P values < 0.07) (Fig. 2B, Supplementary Material, Table S4). The meta-analysis of both cohorts showed P values of 0.005 and 0.012 for unadjusted and adjusted models, respectively. Importantly, among all regressors, ABCD was the only variable that was significantly associated with suicide in the meta-analysis (Supplementary Material, Table S4).
In Cohort 1, the presence of alcohol in blood or urine (Supplementary Material, Table S1) showed no association with ABCD (P = 0.34 by Wilcoxon test). In Cohort 2, ABCD frequencies did not vary among MDD subjects with different severity of depression (Supplementary Material, Table S2) [P = 0.96 by one-way analysis of variance (ANOVA)]. In addition, within each group (all MDD subjects, MDDNon-SU, or MDDSU), there were no associations between ABCD and antidepressant medications (all P values > 0.48 by Wilcoxon test). The latter were determined by positive toxicology for at least one antidepressant at the time of death (Supplementary Material, Table S2). Collectively, these results confirm our previous findings of increased 5-HT2CR editing in the PFC of suicide victims irrespective of the influence of the underlying psychiatric illness (22,39).
5-HT2CR editing shows a robust gene expression signature
We then performed an analysis of association between genome-wide gene expression and ABCD in Cohorts 1 and 2, controlling for age, sex, PMI, brain pH as well as DX in Cohort 2. We first identified genes whose expression correlated with ABCD using a liberal criterion of significance (unadjusted, one-tailed P-value <0.05). Among 10 140 genes that were present in the expression data sets of both Cohorts 1 and 2 (Supplementary Material, Table S5), the expression of 684 genes was positively associated with ABCD in both cohorts [odds ratio (OR) = 2.54, Fisher's exact test overlap P-value <1E−15], and the expression of 346 genes was negatively associated with ABCD in both cohorts (OR = 4.43, Fisher's exact test overlap P-value <1E−15) (Fig. 3, Supplementary Material, Table S6). Collectively 1030 genes were significantly associated with ABCD in both cohorts. We denoted this set of genes as ‘ABCD expression signature’ (or simply ‘ABCD signature’) genes (Supplementary Material, Table S7).
Signature genes are involved in neurotransmission and are linked to mental disorders
We then performed a functional annotation analysis of the ABCD signature genes using MetaCore (www.genego.com), and detected highly significant enrichment for neuron- and nervous system-related terms in all examined ontology databases (Supplementary Material, Table S8). Specifically, the ABCD signature was significantly enriched for genes involved in synaptic transmission and neurotransmitter secretion as well as genes localized within synapse and synaptic vehicle membranes, and genes associated with serious mental illnesses (such as SZ, MDD and BPD). The enrichment for neural-related categories was detected mostly among genes that were positively associated with ABCD, whereas among genes that were negatively associated with ABCD, a number of inflammation and immunity related categories were enriched (Supplementary Material, Table S9).
We detected a strong enrichment of ABCD signature genes for genes involved in both GABA-ergic and glutamatergic synaptic transmission (Supplementary Material, Tables S8 and S9). In particular, the expression of GAD1 (the major GABA-synthesizing enzyme, glutamic acid decarboxylase), SLC32A1 (a vesicular inhibitory acid transporter that is involved in GABA and glycine uptake in synaptic vesicles), GABRA1, GABRA5, GABRB2, GABRB3, GABRD, GABRG2 and GABARAPL1 (multiple subunits of GABAA receptors and GABAA receptor-associated protein-like 1) as well as GRIA4, GRIK2, GRINA and GRM7 (different subunits of ionotropic glutamate receptors and the metabotropic receptor 7) was positively associated with ABCD.
In addition to factors related to specific neurotransmitter systems, synaptic transmission is prominently regulated by presynaptic factors controlling the vesicle-mediated release of neurotransmitters (40). In line with these observations, the ABCD signature included numerous genes that encode various components of synaptic vesicles such as N-ethylmaleimide sensitive fusion protein (NSF), NSF attachment receptor (SNARE) proteins and other vesicle proteins, including synaptosomal-associated protein 25 kDa (SNAP25), synaptobrevins (VAMP1, VAMP2, VAMP5, VAMP8), synapsin (SYN2), syntaxins (STX1B, STX6, STX12), syntaxin binding (or MUNC18) proteins (STXBP1, STXBP2), synaptogyrins (SYNGR1, SYNGR3), synaptotagmins (SYT1, SYT11, SYT13), clathrin proteins (CLTB, CLTC), the Rab family of small GTPases (RAB2A, RAB11A), and double C2-like domain-containing protein alpha (DOC2A).
As mentioned above, 5-HT2CR is coupled to several G proteins that stimulate intracellular responses via effector molecules which influence the activity of ion channels (5). In line with these data, the ABCD signature included genes that encode several G proteins (GNAO1, GNAI1, GNA12, GNB5, GNG3, GNG5, GNG11) as well as multiple (mostly voltage-gated) potassium channels (KCNC1, KCNC2, KCNC4, KCNA1, KCNS3, KCNQ5, KCNF1, KCNK4, KCNJ12, KCTD8, HCN1) and a voltage-gated potassium channel-interacting protein (KCNIP4).
The ABCD signature also included many genes that encode proteins involved in Ca2+ signaling, such as calmodulin (CALM2, CALM3) and Ca2+/calmodulin-dependent protein kinases (CAMKK1, CAMK2D, CAMK2G, MKNK2). Notably, calmodulin has been reported to bind to the proximal region of the 5-HT2CR C-terminus upon receptor activation by 5-HT, which is critical for G protein-independent, arrestin-dependent ERK1,2 signaling mediated by 5-HT2CR (41).
Although in most cases, expression of genes with neuronal functions was found to be positively associated with ABCD, a variety of neuronal genes within the signature showed a negative association (e.g. GNA12, GNG5, GNG11, PLCD1, MKNK2, VAMP5), implying a complex relationship between editing and neuronal functional pathways. Last but not least, the ABCD signature included ADARB1, which encodes ADAR2, but did not include ADAR, which encodes ADAR1 or ADARB2, which encodes ADAR3. ADAR3 is a brain-specific protein whose sequence is closely related to ADAR2, but whose activity has not yet been demonstrated (42). Expression of ADARB1 was positively associated with ABCD.
A recent report identified 41 blood mRNA biomarkers of suicidality using a combination of different designs and approaches (43). This set of blood biomarker genes was significantly enriched for genes related to inflammatory disease and inflammatory response categories. Interestingly, 4 of these 41 genes were present in the ABCD signature (GBP1, ATP6V0E1, EPHX1, TNFSF10) and showed a negative association with ABCD. This is in line with our finding that the set of genes that was negatively associated with ABCD is highly enriched for inflammation and immunity related ontologies.
Distinct cell-type specificity of ABCD signature genes
In our recent study, we compared DNA methylation patterns between Fluorescence-Activated Cell Sorting (FACS)-sorted neuronal and non-neuronal (glial) cells using sister-aliquots of the OFC specimens obtained from six non-suicide subjects from Cohort 1 (44). Classification of the differentially methylated CpG sites into those undermethylated in neurons (neuronal type) and those undermethylated in non-neuronal cells (glial type), combined with findings of others that methylation within control elements typically negatively correlates with gene expression, yielded large sets of predicted neuron-specific and glia-specific genes (5098 and 6821 genes, respectively). These sets of genes were in excellent agreement with the available direct measurements of gene expression in human and mouse. In order to explore cell-specificity of the ABCD signature genes, we investigated the overlap between the signature genes and these DNA methylation data. We detected a highly significant overlap between those signature genes that were positively associated with ABCD (684 genes) and the neuron-specific gene set (OR = 1.988, P-value < 1E−15 by Fisher's exact test), whereas there was no significant overlap between these signature genes and the glial-specific genes. In contrast, the group of the signature genes that were negatively associated with ABCD (346 genes) showed a non-random overlap with the glial-specific gene set of (OR = 3.312, P-value < 1E−15) but not the neuronal-specific set (Supplementary Material, Table S10).
We next compared the ABCD signature genes with the cell type-specific gene-expression patterns derived from the analysis of the human transcriptome. We applied an in silico tissue dissection statistical approach that is based on ‘weighted gene co-expression network analysis’ (WGCNA) (45). The WGCNA identifies biologically relevant patterns in expression datasets by grouping genes into modules with strongly covarying patterns across the samples. Consequently, WGCNA can distinguish gene-expression networks associated with specific cell types (e.g. neurons or glia) that are present in a heterogeneous sample (e.g. whole human cortex) thanks to the distinct transcriptional profiles of these cell types and variation in the relative proportions of neurons and glial cells across samples (46). In particular, WGCNA co-expression modules can be examined for cell-type specificity using cell type-enriched gene sets identified in previous studies (47–53). We have previously applied WGCNA to the postnatal sub-cohort (N = 220) of the gene expression microarray dataset from the BrainCloud study which analyzed PFC samples obtained from a large collection of autopsy brain specimens (Gene Expression Omnibus; accession GSE30272) (54). Consistent with previous studies in the brain tissues (46,51), the WGCNA yielded well-defined gene co-expression modules (N = 77) with specific anatomical distributions frequently related to specific cell types (44). We detected numerous modules consisting of genes with enhanced expression in neurons or in different types of glia (microglia, astrocytes and oligodendrocytes) (44). A comparison of the coexpression modules with the ABCD geneexpression signature showed that the signature genes positively associated with ABCD (N = 684) were mostly over-represented in the modules related to glutamatergic or GABAegric (e.g. parvalbumin-expressing) neurons (adjusted P-values <0.01), whereas those signature genes that were negatively associated with ABCD (346 genes) were over-represented in microglia and astrocyte modules (all adjusted P-values <0.05) (Supplementary Material, Table S11). We did not detect an enrichment of the negatively associated signature genes in modules related to oligodendrocytes. Instead and in line with our gene ontology findings which linked the negatively associated signature genes with categories related to inflammation, 5 of the 8 modules that were enriched for the negatively associated signature genes were related to microglia.
Signature genes are over-represented among genes which show association between expression and DNA methylation among different individuals
Animal studies indicate that 5-HT2CR editing (and therefore, function) can be regulated by environmental factors (4). Because environment can influence gene expression via epigenetic processes, we asked if the ABCD signature was enriched for genes whose expression is linked to epigenetic modifications. Although there are many epigenetic marks that potentially could be studied, DNA methylation is the major mark that has already been demonstrated to be stable postmortem (55). Thus, we investigated the relationship between DNA methylation and expression of the ABCD signature genes.
The present results indicate that the majority of the ABCD signature genes are expressed in neurons. In addition, a dramatic difference in the DNA methylation landscapes between neuronal and glia cells has been recently reported (44,56). Therefore, we assessed DNA methylation within the PFC neuronal population. Neuronal nuclei were isolated from the sister aliquots of the OFC specimens from Cohort 1 (N = 51), which were used in our 5-HT2CR editing and RNA-Seq studies described above. We employed Illumina HM450K array technology, which examines the methylation status of ∼485 000 individual cytosine positions (mostly CpG sites). Although the array contains only ∼2% of the total CpG sites present in the genome, it spans the entire genome covering 99% of the RefSeq protein-coding genes with multiple probes per gene, 96% of the CpG islands from the UCSC database, high and low CpG density promoters and predicted enhancers.
We first identified those CpG sites in the HM450K array that belonged to the 10 140 genes whose expression were analyzed in both cohorts and were used to delineate the ABCD signature (see Supplementary Material, Supplementary Methods for CpG site-to-gene annotation). Next, regression analysis was applied to determine which of these 10 140 genes had at least one CpG site whose methylation level correlated with the expression of the respective gene in Cohort 1. To avoid spurious associations, we restricted the analysis to those sites which showed at least 10% range of variation in the beta values (that is, the sites for which the difference between the maximum and minimum beta value was at least 0.1). We detected a set of 551 such genes (hereinafter, methylation-linked genes). Among the methylation-linked genes, 63 were present in the ABCD signature, and 56 of these 63 genes were positively associated with ABCD (Supplementary Material, Table S12). Thus, although there was no significant overlap between the signature as a whole (1080 genes) and the set of the methylation-linked genes (OR = 1.150, P = 0.173), there was a significant overlap between the 684 positively associated signature genes and this set (OR = 1.644, P = 8.584E−4). No significant overlap was detected with negatively associated signature genes (OR = 0.340, P = 0.999). Because DNA methylation was measured in neurons, this observation is compatible with our finding that the positively associated signature genes are preferentially expressed in neurons. Thus, the ABCD signature (or at least its neuronal component) is enriched for genes whose expression is linked to DNA methylation.
We also tested the overlap between the ABCD signature and a list of genes involved in epigenetic regulation (57–60), whose expression was detected in both cohorts. Of the 69 such genes, only 3 were present in the signature: HDAC1, SMYD2 and SMYD3 (Supplementary Material, Table S13). This overlap, however, was not significant (OR = 0.40, Fisher's exact test P-value = 0.98).
Association between editing and gene expression is altered in suicide
We then asked whether or not the relationship between editing and gene expression differs in people who died of suicide compared with individuals who died of non-suicide causes. We first developed an overall ABCD expression signature score for each individual. The score was computed as the difference between the average standardized expression of the signature genes positively associated with ABCD and the average standardized expression of the signature genes negatively associated with ABCD (see Materials and Methods). Accordingly, for each individual, the signature score summarized the expression pattern of all signature genes, taking into account the sign of their association with ABCD. A high signature score for a specific individual indicated high expression of the signature genes positively associated with ABCD and low expression of the signature genes negatively associated with ABCD. Conversely, a low signature score indicated low expression of the genes positively associated with ABCD and high expression of the genes negatively associated with ABCD.
We then compared the signature scores between the suicides and non-suicides in each cohort. As expected, we observed significant association between the signature score and ABCD in non-suicide groups in both cohorts (all P-values < 1E−4) (Fig. 4). In a striking contrast, among the suicide victims from both cohorts, there was no association between the score and ABCD (P-values > 0.462). In other words, the range of variation for both ABCD frequencies and signature scores was dramatically lower in the suicide compared with the non-suicide groups.
We then attempted to identify the individual genes that were associated with ABCD in non-suicide but not in suicide subjects. To this end, we performed an association analysis between ABCD and gene expression in both cohorts introducing ‘Suicide’ as an interaction term. This analysis yielded no significant overlap between the cohorts, probably because much larger samples are needed to detect interactions (61). However, because the signature score captures the expression patterns for all ABCD signature genes, the differences in the patterns of association between the suicide and non-suicide groups for the majority of the genes in the signature were very similar to each other and to the overall pattern (Fig. 5 and Supplementary Material, Fig. S3).
RNA editing is one of the major mechanisms by which environmental signals overwrite genetically encoded information to modify gene function and regulation, particularly in the brain, thus providing a possible interface between environmental inputs and behavior. Higher overall editing levels have been reported for the human brain compared with the brains of other mammals including non-human primates (8,11), suggesting that RNA editing might underpin higher-order brain function and complex behaviors, and could be involved in pathophysiology of psychiatric and neurological diseases.
Previous studies have shown that 5-HT2CR editing is altered in the PFC of suicide victims (22,23,29–31), and we confirmed these findings in the present work. Hence, we speculate that A-to-I editing in general, and 5-HT2CR editing in particular, could act as a homeostatic mechanism that dampens the effect of a wide range of environmental and genetic factors that impair neuronal function. This hypothesis implies that 5-HT2CR editing is linked to the signaling pathways that are important for neuronal activity and consequently to the expression of genes whose products are involved in these pathways. Similar hypotheses have been recently proposed on the basis of a comparative study of A-to-I editing in autistic and non-autistic cerebella (62). We further speculated that the relationship between 5-HT2CR editing and gene expression would differ between individuals who died of suicide versus those who died of non-suicide-related causes.
The present study aimed to test these hypotheses in two cohorts of individuals who died of suicide or other causes. Although the cohorts were, out of necessity, relatively small, the results show that 5-HT2CR editing is associated with a robust gene expression signature, i.e. a large set of genes whose expression was significantly correlated with editing in both cohorts. This signature set was significantly enriched for genes involved in synaptic function and genes associated with serious mental illnesses, such as SZ, MDD and BPD, as well as suicide (Supplementary Material, Tables S8 and S9). Psychological autopsy studies show that >90% of all individuals who die by suicide meet criteria for at least one psychiatric disorder, most commonly MDD, BPD, SZ or substance use disorder (32).
Our analysis also demonstrated that the signature set was enriched for genes that are preferentially expressed in neurons and for those neuronal genes whose expression is negatively correlated with their DNA methylation. In the brain, DNA methylation within promoters and distal regulatory regions is associated with repression of gene expression (44,56). In addition, although DNA methylation marks in postmitotic cells may persist throughout the cellular lifetime, the methylation status is a flexible substrate for epigenetic modification that can be altered by cellular response to environmental stimuli (63). Thus, our findings imply a link between environmental impact and 5-HT2CR editing through chromatin remodeling and the associated gene expression changes.
The key finding of this study is that the association between editing and gene expression is disrupted in suicide, suggesting a dysregulation of the putative homeostatic function of 5-HT2CR editing in the PFC of suicide victims. Specifically, compared with the non-psychiatric controls and especially to the MDD without suicide group, the suicide group included no individuals with low levels of editing and low ABCD signature scores (Figs 4 and 5). We speculate that in the brain of suicide victims, a mechanism that enhances the function of 5-HT2CR through lowering the level of editing and hence prevents neuron malfunction is impaired. Elucidation of the nature of this putative mechanism(s) is a goal for subsequent work.
The editing-associated signature included ADARB1, which encodes ADAR2, but not ADAR which encodes ADAR1. A recent study has reported that the dysfunctional isoform of ADAR2 is over-represented in postmortem cerebella from individuals with autism (62). Thus, the pathways related to the regulation of ADAR2 might be important for gene-expression mediated crosstalk between neuronal function and 5-HT2CR editing. These, still virtually uncharacterized ADAR2-related pathways (64) could be compromised in suicide. Indeed, as shown here, the apparent dysregulation of editing in the PFC of suicide victims is characterized by a notable increase in the frequency of the ABCD variant the formation of which requires editing by both ADAR1 and ADAR2 because site D is primarily edited by ADAR2 (65).
The most obvious shortcoming of the present work appears to be that the majority of the MDD subjects in Cohort 2 were receiving antidepressant medications. Therefore it was difficult to ascertain whether medication, rather than an innate pathophysiology, accounted for the differences in editing between MDDNon-SU and MDDSU. Multiple lines of evidence, however, argue against the possibility that suicide-related differences in editing observed in our study could be caused by medication: (1) Similar to the results obtained in earlier animal studies (reviewed in 4), a recent report in mice that employed MPS (and therefore, achieved high precision in assessing 5-HT2CR editing levels) showed that chronic treatments with antidepressant drugs resulted in an increase in 5-HT2CR mRNA editing (66). Specifically, in that study increased editing at the A and B sites was observed in the striatum and hippocampus following fluoxetine and only in the hippocampus following amitriptyline treatment. Neither drug influenced editing in the neocortex; in addition, no editing alterations were observed in any brain region following treatment with antipsychotic medications (olanzapine or clozapine); (2) In the present study, within each group in Cohort 2 (all MDD subjects, MDDNon-SU, or MDDSU), there was no association between ABCD and antidepressant medications. Moreover, the percentage of subjects that showed positive toxicology for at least one antidepressant medication at the time of death was higher in the MDDNoSuic (75%) than in the MDDSuic (60%) group; however, more editing was detected in the MDDSuic than in the MDDNoSuic; (3) In our previous work on 5-HT2CR editing in SZ and BPD patients with and without suicide, the history of antidepressant treatment among suicide victims had no discernible impact on editing (22). Also, no differences were observed between non-psychiatric controls and SZ or BPD patients that did not commit suicide, although more than half of the non-suicide BPD patients received antidepressant medication, and all non-suicide SZ patients received antipsychotic medications; (4) In the present study, we observed difference in editing between suicide victims and subjects without any antemortem DX in Cohort 1. Although 6 of the 22 suicide subjects in this cohort had a history of depression, only one suicide subject showed positive toxicology for medications (olanzapine and clozapine). Based on this collective evidence, we believe that the editing alterations in suicide subjects that we identified in this work are unlikely to be related to their treatment with antidepressant or antipsychotic medications.
To summarize, the results of this work imply that the postulated homeostatic function of 5-HT2CR editing is dysregulated in the brains of suicide victims. The resulting increased editing (or more precisely, relative paucity of unedited variants) could be a synergistic effect of suicide-associated mutations and environmental stimuli. The relatively small number of samples analyzed in this report does not allow us to examine this possibility. In studies of suicide, which is a rare and unpredictable event, obtaining postmortem samples is an obvious challenge. We attempted to partially overcome this problem by studying the link between editing, gene expression and DNA methylation. Further, more detailed and accurate analysis of the multiple connections between genetic and molecular phenotypic signatures of suicide seems to have potential to unravel the molecular mechanisms underpinning this complex behavior and ultimately develop effective interventions.
MATERIALS AND METHODS
Subjects and brain regions
Human postmortem brain specimens from the PFC of Caucasian individuals were provided by the Brain Collection of Dr Yasmin Hurd (Icahn School of Medicine at Mount Sinai, New York, NY) (Cohort 1) (67,68) and by the University of Pittsburgh Brain Tissue Donation Program (Pittsburgh, PA) (Cohort 2) (69,70). Detailed description of the cohorts and brain regions is provided in Supplementary Methods and Supplementary Material, Tables S1 and S2. All brains were analyzed for adequate brain pH [6.8 ± 0.2 and 6.6 ± 0.3, mean ± standard deviation (SD) for Cohorts 1 and 2, respectively].
In Cohort 1, 6 of the 22 suicide subjects had a history of depression. The NC subjects did not have any antemortem DX. At the time of death, all subjects showed negative toxicology for opiates and other illicit drugs. Three control subjects and five suicide subjects were positive for alcohol in blood or urine, and therapeutic agents were detected in one control and one suicide subject.
In Cohort 2, control subjects showed negative toxicology for illicit drugs and antidepressant drugs, while 18 MDDNon-SU and 6 MDDSU individuals showed positive toxicology for at least one antidepressant medication at the time of death.
RNA extraction, quality control and analysis of 5-HT2CR RNA editing
Total RNA was extracted using Invitrogen ToTally RNA Kit (Cohort 1) or TRIzol reagent (Cohort 2) (both reagents from Life Technologies). The quality of RNA was assessed using Agilent Bioanalyzer. RNA integrity number (RIN) for all specimens was >6.4 (7.6 ± 0.6 and 8.1 ± 0.6, mean ± SD for Cohorts 1 and 2, respectively). Complimentary DNA (cDNA) was synthesized from equal quantities of RNA from each subject (1 μg of RNA per 10 μl of reaction) using High Capacity cDNA Reverse Transcription kit (Life Technologies).
5-HT2CR RNA editing was assessed using MPS of PCR-amplified fragments in the region of editing as described in refs (39,71) and detailed in Supplementary Methods. Because inosine is read as guanine during reverse transcription (6), the sequencing enables to assess the pattern and the level of editing in each experimental sample. Specifically, MPS provides the number of reads for each of the 32 5-HT2CR mRNA editing variants in each subject. The relative frequency of each variant in each subject was calculated as the ratio of the number of reads detected for a particular variant to the total number of reads for all variants in the subject.
Transcriptional profiling of Cohort 1 was performed using mRNA-Seq. For each subject, the same RNA preparations were used for both RNA editing and RNA-Seq assays. Libraries were prepared from 1 μg of total RNA using Illumina's TruSeq RNA library preparation kit V2, and sequencing was performed using an Illumina HiSeq2000. After de-multiplexing of each library pool according to its barcode sequence, from 8 to 20 million (M) (mean 15.1 M) reads per sample were available for analysis. Alignment, assembly and quantification of the data were performed using the TopHat 2.0.8b and Cufflinks 2.1.1 (72). Only transcripts with the maximum value of the fragments per kilobase of transcript per million mapped reads (FPKM) >2 among all subjects were analyzed in the study. Gene expression profiles of Cohort 2 specimens have been previously assessed using Illumina microarray chip Human v3, which contains 48 803 RefSeq probes (73).
For each gene that was represented by several transcripts (Cohort 1) or probes (Cohort 2), only one transcript/probe that resulted in a signal with the highest SD was retained. After this filtering, the data for 13 101 and 24 022 unique gene symbols were retained for Cohorts 1 and 2, respectively. Among these, there were 10 140 gene symbols that were present in both datasets (cohorts), and those genes were used in subsequent association analyses with 5-HT2CR editing and DNA methylation data (see below). The list of genes that were present in both datasets (cohorts) with associated P-values is reported in Supplementary Material, Table S5.
Analysis of association between editing and gene expression
The association between gene expression and editing was assessed for Cohorts 1 and 2 by fitting linear regression models of gene expression level versus editing level. The editing level was represented by the frequency of the ABCD mRNA variant in each individual. The model also included other factors which could affect gene expression as covariates, such as age, sex, PMI, brain pH and DX (suicide or NC in Cohort 1; MDDSU, MDDNon-SU or NC in Cohort 2). The categorical variables (DX and sex) were expanded into k− 1 indicator variables, k being the number of levels observed for that variable.
Calculation of ABCD gene expression signature score
The ABCD gene expression signature score was derived in order to summarize the expression pattern of the ABCD signature genes in each individual. The score allowed us to analyze the relationship between the ABCD signature and the clinical covariates (DX and suicide status). To calculate the score, the expression values for signature genes were first standardized to a mean of ‘0’ and a SD of ‘1’. In Cohort 1, the FPKM expression values (incremented by 1) were log-2 transformed prior to the standardization. The score was then computed for each individual as the difference between the average standardized expression of the signature genes that were positively associated with ABCD and the average standardized expression of the signature genes that were negatively associated with ABCD.
DNA methylation analysis of neuronal cells
The DNA methylation of neuronal cells was measured in Cohort 1 only. Neuronal nuclei were separated from non-neuronal (glial) nuclei using FACS with DNA-binding dye 7-AAD and anti-NeuN neuron-specific antibodies mostly as described in ref. (44) and as detailed in Supplementary Methods.
Genomic scale single-site resolution DNA methylation profiling of neuronal cells was performed using Illumina's Infinium Human Methylation450K array (HM450K) (Illumina Inc.) as described in refs (74,75) and detailed in Supplementary Methods. Methylation level at each individual CpG site in each specimen was measured as ‘beta value’ (β), which is defined as the ratio of the methylated probe intensity to the sum of the methylated and unmethylated probe intensities. The resulting β can range from 0 (no methylation) to 1 (complete methylation). The CpG site-to-gene annotation was performed as described in Supplementary Methods.
Analysis of association between DNA methylation and gene expression
The association between gene expression and DNA methylation was analyzed in Cohort 1 (N = 51) using a linear regression model. Only those 10 140 genes that were detected in both cohorts were analyzed. Age, sex, PMI and brain pH were included in the statistical model. The analysis was performed for cis associations only. Specifically, expression of each gene was correlated with the methylation level at the CpG sites that were annotated for the given gene.
Analysis of difference in editing among disease groups
The differences in the average level of the ABCD mRNA variant among disease groups were tested using the two-sided Wilcoxon test in Cohort 1 and the one-way ANOVA followed by the Tukey's HSD post hoc pairwise comparisons in Cohort 2. Prior to each test, within each Cohort–disease group combination, observations below the 25th percentile minus 1.5 times the interquartile range, and above the 75th percentile plus 1.5 times the interquartile range were identified as outliers (resulting in four outliers in Cohort 1 and four outliers in Cohort 2) and were removed from the analysis.
Suicide propensity modeling
Univariate (one regressor at a time) and multivariate (all regressors: ABCD, age, sex, PMI as well as DX for Cohort 2) logistic regression models were fitted for the probability of suicide against the available regressors in Cohorts 1 and 2.
Linear and logistic modeling, and PCA as well as visualization were performed using R, version 3.0.0 (http://www.R-project.org/). The PCA of the frequencies of 5-HT2CR mRNA editing variants was performed using the ‘prcomp’ function in R. Linear modeling of gene expression data was performed using Limma, version 3.12.1.
This work was funded by VISN3 Mental Illness Research and Education Clinical Center (MIRECC), VA Merit Review Award (BX001829 to S.D.); National Institute of Health (R21MH090352 and R21DA031557 to S.D.; and R01DA015446 to Y.H.); and Hope for Depression Research Foundation grant to S.D. E.V.K. is supported by intramural funds of the US Department of Health and Human Services to National Library of Medicine.
The authors wish to thank Alisa Timashpolsky for superb technical assistance. The work was supported with resources and the use of facilities at the James J Peters VA Medical Center, Bronx, New York.
Conflict of Interest statement. Authors have no competing financial interests in relation to the work described.