Transcriptome-Wide Identification of G-to-A RNA Editing in Chronic Social Defeat Stress Mouse Models

Emerging evidence suggests that RNA editing is associated with stress, neurological diseases, and psychiatric disorders. However, the role of G-to-A RNA editing in chronic social defeat stress (CSDS) remains unclear. We herein identified G-to-A RNA editing and its changes in the ventral tegmental area (VTA), a key region of the brain reward system, in CSDS mouse models under emotional stress (ES) and physiological stress (PS) conditions. Our results revealed 3812 high-confidence G-to-A editing events. Among them, 56 events were significantly downregulated while 23 significantly upregulated in CSDS compared to controls. Moreover, divergent editing patterns were observed between CSDS mice under ES and PS conditions, with 42 and 21 events significantly upregulated in PS and ES, respectively. Interestingly, differential RNA editing was enriched in genes with multiple editing events. Genes differentially edited in CSDS included those genetically associated with mental or neurodevelopmental disorders, especially mood disorders, such as FAT atypical cadherin 1 and solute carrier family 6 member 1. Notably, changes of G-to-A RNA editing were also implicated in ionotropic glutamate receptors, a group of well-known targets of adenosine-to-inosine RNA editing. Such results demonstrate dynamic G-to-A RNA editing changes in the brain of CSDS mouse models, underlining its role as a potential molecular mechanism of CSDS and stress-related diseases.

Introduction

RNA editing is the change of RNA nucleotide sequence that occurs at the transcriptional level (Tan et al., 2017). Canonical RNA editing in mammals includes adenosine to inosine (A-to-I) editing mediated by the adenosine deaminase acting on RNA (ADAR) protein family and cytidine to uridine (C-to-U) mediated by the Apolipoprotein B mRNA Editing Catalytic Polypeptide-like (APOBEC) protein family (Gu et al., 2012). RNA editing has been found in various organisms, including animals (Tan et al., 2017), plants (Small et al., 2020), and humans (Christofi and Zaravinos, 2019), and is involved in a variety of physiological and pathological processes. Recent studies have suggested that RNA editing significantly changes in response to stress, neurological diseases and psychiatric disorders (Kawahara et al., 2004; Breen et al., 2019; Dick et al., 2019), suggesting that it may be involved in the molecular mechanisms of these pathological processes.

Non-canonical RNA editing refers to change of RNA sequence other than A-to-I or C-to-U canonical RNA editing. Compared to canonical RNA editing, non-canonical RNA editing remains relatively poorly understood. G-to-A RNA editing has been reported in the heterogeneous nuclear ribonucleoprotein K (HNRNPK) gene in colorectal cancer and the glutamate ionotropic receptor kainate type subunit 3, GRIK3) gene in the fetal human brain (Nutt et al., 1994, 7; Klimek-Tomczak et al., 2006). G-to-A editing in WT1 transcripts (Niavarani et al., 2015) was reported to be mediated by the deaminase activity of APOBEC3A. Nevertheless, the role of G-to-A editing in physiological and pathological processes is largely unclear.

Social defeat stress models are important tools to study post-traumatic stress disorder (PTSD), depression, and other stress-related diseases. Recent studies have shown that social defeat stress could induce changes in RNA editing (Dick et al., 2019). Dick et al. (2019) demonstrated dynamic regulation of A-to-I editing in the prefrontal cortex and basolateral amygdala after chronic social defeat (CSDS) in mice. However, the role of non-canonical G-to-A RNA editing in CSDS remains to be elucidated.

Chronic social defeat stress could be induced in mice under emotional stress (ES) condition by witnessing traumatic events and insulating effects of physical stress (PS) (Warren et al., 2013). Interestingly, both similarities and differences in brain transcript expression profiles were found between CSDS mice under ES and PS conditions in Warren’s study. Therefore, G-to-A RNA editing in CSDS mouse models under different stress conditions is still to be studied. The current study herein identified G-to-A RNA editing from RNA sequencing (RNA-Seq) data of the mouse ventral tegmental area (VTA) at the transcriptome level. Our findings revealed that dynamic changes of G-to-A RNA editing were associated with CSDS, and could also contribute to the difference between CSDS models under ES and PS conditions.

Materials and Methods

Collection of RNA-Seq Data

Raw RNA-Seq read data were downloaded from NCBI Gene Expression Omnibus (GEO) database¹. The samples used for RNA editing event discovery consisted of adult C57BL/6J male mouse brain VTA from CSDS mouse models under ES or PS conditions, and controls (N = 3 per group) (GSE36005) (Warren et al., 2013). As described Warren’s study, naïve C57BL/6J mice were assigned to either ES or PS conditions. ES mice were exposed to witnessing the social defeat of a PS mouse attacked by an aggressor CD-1 mouse from a safe adjacent compartment. Another set of VTA samples from 44 adult C57BL/6 mice were analyzed independently to generate a reference dataset, and used to help identify high-confidence RNA editing events (GSE89692) (Peña et al., 2017, 2019).

Processing of RNA-Seq Data

The raw sequencing data obtained above were first analyzed by FASTQC for quality control. Reads that passed quality control were aligned and mapped to the mouse genome (UCSC mm10) using RNA STAR (version 2.7.0e) (Dobin and Gingeras, 2015), and generating alignment files in BAM format. SamTools (version 1.9) was used to filter the reads by removing optic duplications in the BAM files (Li et al., 2009), and only reads uniquely mapped to the mouse genome were kept. Base quality score recalibration was then performed with the resulting BAM files by using GATK (version 4.1.3) and following the best practices workflows recommended by the documentation (Van der Auwera et al., 2013).

Identification of High-Confidence G-to-A RNA Editing Events

RNA-Seq data from CSDS and control mice (GSE36005) were used for discovery of RNA editing events. Single nucleotide variants (SNV) were called from the GATK re-calibrated BAM files by using VarScan (version 2.4.3) (Koboldt et al., 2012). The variant calling criteria were set as follows: base quality ≥25, total sequencing depth ≥10, alternative allele depth ≥2 and alternative allele frequency (AAF) ≥1%, and possible false positive SNVs were filtered and removed by using the fpfilter command implemented by VarScan with default parameters. SNVs annotated using the Ensembl Variant Effect Predictor (VEP)² (McLaren et al., 2016). SNVs that met any of the following conditions were further filtered and removed unless annotated as known RNA editing events in the REDIportal V2.0 database³ (Mansi et al., 2021): (1) located in simple repeats or homopolymer runs ≥5 nucleotides (nt); (2) located in the mitochondria; (3) located within 6 nt from splice junctions; (4) located within 1 nt from insertion-deletions (INDELs); (5) annotated as known variants in the dbSNP database; or (6) more than 90% of the samples of controls and CSDS had an AAF equal to 100% or between 40 and 60%. RNA-Seq data from another 44 adult mice (GSE89692) were used as an additional dataset to generate a list of reference G-to-A RNA editing events, and help identify high-confidence ones. The same analysis workflow used in the RNA editing event discovery was applied. Only G-to-A SNVs that were located in genic regions, and were detected in at least 2 samples or replicated in the additional set of 44 mouse VTA samples with editing levels ≥1% were kept, which were defined as high-confidence G-to-A RNA editing events.

Quantification of Gene Expression

FeatureCounts was used to obtain the pseudo-counts of gene expression from the BAM files (Liao et al., 2014), and transcripts read per thousand bases per million mapping (TPM) were then calculated for each gene using EdgeR (version 3.7) (McCarthy et al., 2012).

Principal Component Analysis

Principal component analysis (PCA) of G-to-A RNA editing events was performed using the function Prcomp in R (version 3.6.3), and visualized using the ggplot2 (version 2.2.1) package.

Enrichment Analysis of Gene Functions and Pathways

In order to understand the potential biological effects of RNA editing, the analysis of the Gene list enrichment analysis of gene ontology (GO) and kyoko encyclopedia of genes and genomes (KEGG) pathways were conducted using Enrichr⁴ (Kuleshov et al., 2016). Enrichment of with by FDR < 0.05 as the cutoff.

Statistical Analysis

Levels of RNA editing or gene expression were compared using Student’s t-test between groups. P-values < 0.05 were considered to be significant. The correlation coefficient (r) between RNA editing levels was calculated using Spearman correlation analysis. Fisher’s exact test was used to analyze enrichment of differentially edited genes.

Results

G-to-A RNA Editing Events Identified in Mouse VTA Transcriptome

Our bioinformatic analysis identified a total of 3812 high-confidence G-to-A RNA editing events in 2553 genes in the mouse VTA RNA-Seq data (Figure 1 and Supplementary Table 1). These events, with editing levels of these events ranging from 1 to 100%, were widely observed across all chromosomes. These identified high-confidence RNA editing events consisted of various functional categories, including 1492 missense, 1333 3′-untranslated region (UTR), 726 synonymous, 121 5′-UTR, 76 non-coding transcript exons variants, 63 stop-gain, and 1 stop-lost variants (Figure 1B). SIFT predicted 754 (50.5%) of the missense events could have potential impact on the protein function or structure (Figure 1C).

FIGURE 1

Figure 1. G-to-A RNA editing events identified from adult mouse VTA transcriptome. The red dots denotes the events observed across all chromosomes in the mouse genome (A). The blue dots show mean expression levels of genes. And the lines denote interaction between the G-to-A RNA editing events, with the blue-to-red gradient indicated the correlation co-efficient r. The interaction among the top 100 frequently observed editing events is shown. (B) The G-to-A RNA editing events result in various types of mRNA variants. (C) About half of these missense events are predicted by SIFT to possibly be deleterious on the encoded proteins.

CSD-Associated G-to-A RNA Editing Events

As shown in Figure 2, comparison of G-to-A RNA editing between CSDS and controls revealed that, 2729 (71.6%) were shared by both CSDS and controls, 865 (20.8%) were observed only in CSDS but not in controls, and 871 (24.2%) only in controls but not in CSDS. T-test identified a total of 79 G-to-A RNA editing events with significantly different editing levels, among which 56 were downregulated and 23 were upregulated in CSDS compared to controls, respectively (Figure 3A). These CSDS-associated events consisted of 36 in 3′-UTR, 2 in 5′-UTR, 4 in non-coding transcript exons, and 20 in coding regions including 20 missense and 17 synonymous variants, with no significant expression changes in most of the genes with differential editing (Supplementary Table 2). 22 of the genes with differential G-to-A editing in CSDS had been reported to be associated with psychiatric diseases, and another five of them were associated with other neurological diseases. PCA with these 79 differential editing events showed that CSDS samples clustered separately from controls, with 55.62% contribution from PC1 to the total variance (Figure 3B).

FIGURE 2

Figure 2. Comparison of genomic locations of G-to-A RNA editing events in controls and CSDS.

FIGURE 3

Figure 3. Differential G-to-A RNA editing events between CSDS and controls. (A) A total of 79 events are differentially edited between CSDS and controls, including 56 and 23 are upregulated and downregulated, respectively. (B) Principle component analysis of the 79 events differentially edited between CSDS and controls.

The three G-to-A RNA editing events that were the most significantly downregulated in CSDS compared to controls were observed in FAT atypical cadherin 1 (Fat1:chr8:45026388, P = 8.7 × 10^–10), TMF1 regulated nuclear protein 1 (Trnp1:chr4:133497722, P = 1.7 × 10^–5), and prothymosin alpha (Ptma:chr1:86529182, P = 3.4 × 10^–5) (Figure 4A). Notably, the Fat1 RNA editing at Fat1:chr8:45026388, which produced a missense change p.V280M (Gtg- > Atg), and another Fat1 missense event at Fat1:chr8:45025105 (aGc/aAc, p.S2373N) were observed in controls only. Both of the two Fat1 G-to-A changes were predicted to be deleterious to the protein function (Supplementary Table 2). In contrast, the top three G-to-A RNA editing events upregulated in CSDS were found in solute carrier family 6 member 1 (Slc6a1:chr6:114316424, P = 0.0036), tripartite Motif Containing 2 (Trim2:chr3:84190826, P = 0.008), and apolipoprotein D (Apod:chr16:31297353, P = 0.0082) (Figure 4B). In addition, our results also implicate that glutamate ionotropic receptor NMDA type subunit 2b (Grin2b) and subunit associated protein 1 (Grina) could be the two most evident glutamate receptors subjected to G-to-A RNA editing in CSDS (Grin2b:chr6:135720485, P = 0.0178; and Grina:chr15:76248349, P = 0.0287, respectively) (Supplementary Table 2). These RNA editing events, either differentially edited or exclusively detected in either controls or CSDS, suggested that changes of G-to-A RNA editing could be associated with CSDS.

FIGURE 4

Figure 4. The events that are the most significantly associated with CSDS. The top three upregulated events (A), and top three downregulated events (B) in CSDS compared to controls are shown.

Functional Enrichment in CSDS-Associated G-to-A RNA Editing

To help understand the biological impact of G-to-A RNA editing in CSDS, enrichment analysis of genes with RNA editing enriched in either controls or CSDS was further compared. Intriguingly, the results in Figure 5 showed stronger functional enrichment in CSDS than in controls. GO analysis showed that biological processes enriched in CSDS included myelination, negative regulation of macroautophagy, cellular protein modification, Ras protein signal transduction, Golgi vesicle transport, small GTPase mediated signal transduction, and protein localization to plasma membrane (Figure 5A). The molecular functions enriched in CSDS were mainly related to cadherin binding, ubiquitin-related enzyme activity, sodium channel activity, Rab GTPase binding, and RNA binding (Figure 5B). The cellular components enriched in CSDS were mainly related to ficolin-1-rich granule, cullin-RING ubiquitin ligase complex, Golgi subcompartment and neuron-specific components such as main axon (Figure 5C). The KEGG pathways enriched in CSDS were mainly involved in insulin signaling, cAMP signaling, endocytosis, ubiquitin mediated proteolysis, glutamatergic and GABAergic synapse, regulation of actin cytoskeleton as well as neurological diseases such as nicotine addiction (Figure 5D).

FIGURE 5

Figure 5. Gene ontology and KEGG pathways erniched in gene with CSDS-associated G-to-A RNA editing. The items with the most significant P-values are shown for (A) biological processes, (B) molecular functions, and (C) cellular components, as well as (D) KEGG pathways enriched in either CSDS or controls. CSDS, chronic social defeat stress.

Divergent G-to-A RNA Editing Patterns Between Different CSDS Models

Although both ES and PS mice were reported to exhibit similar defects of depression- and anxiety-like behavior, ES mice experienced CSDS by only witnessing emotional events without direct physical interaction. Previously Warren’s study indicated possible differences in transcript expression between the two CSDS models. We thus went on to compare the G-to-A RNA editing patterns between the two CSDS paradigms.

Figure 6 showed that among these RNA editing events, 1976 (55.0%) were shared by both ES and PS, 787 (20.8%) were observed only in PS but not in ES, and 871 (24.2%) only in ES but not in PS. Furthermore, Student’s t-test identified 63 events with significantly different editing levels between ES and PS, among which 42 were upregulated in ES and 21 were upregulated in PS (Figure 7A and Supplementary Table 3). Five of the genes differentially edited between ES and PS were associated with psychiatric diseases, and another six of them were associated with other neurological diseases. As shown in Figure 8, the top three events that significantly upregulated in ES were found in cullin 4A (Cul4a:chr8:13133758, P = 8.09 × 10^–9), c-Maf inducing protein (Cmip:chr8:117436681, P = 1.01 × 10^–7), and multivesicular body subunit 12b (Mvb12b:chr2:33730066, P = 9.24 × 10^–7, respectively). The top three events that significantly upregulated in PS were found in SSX family member 2 interacting protein (Ssx2ip:chr3:146428062, P = 0.001), DAZ interacting zinc finger protein 3 (Dzip3: chr16:48924538, P = 0.001), and dystrophin related protein 2 (Drp2:chrX:134454263, P = 0.001). PCA with these 63 events showed that ES cluster separately from PS, with higher intra-group variance in ES, with PC1 contributing 80.2% to the total variance, (Figure 7B). These events, either enriched in ES or PS, indicated that divergent patterns of G-to-A RNA editing might contribute to the underlying different mechanisms in the two CSDS models.

FIGURE 6

Figure 6. Comparison of genomic locations of G-to-A RNA editing events between ES and PS.

FIGURE 7

Figure 7. Differential G-to-A RNA editing events between PS and ES. (A) 63 events showed differential editing levels between CSDS and controls, including 42 and 21 are upregulated in PS and ES, respectively. (B) Principle component analysis of these Differential G-to-A RNA editing between PS and ES.

FIGURE 8

Figure 8. The top G-to-A RNA editing with significant difference between PS and ES. The top three upregulated events in ES (A), and top three downregulated events (B) in PS are shown.

Functional Divergence of G-to-A RNA Editing Enriched in ES or PS

Enrichment analyses of gene functions and pathways were further evaluated in genes with G-to-A RNA editing enriched in either ES or PS (Figure 9). The results demonstrated that a number of these gene functions were found in both CSDS models. Biological processes such as neuron projection development, negative regulation of phosphatase activity, and protein localization to membrane, molecular functions such as RNA binding and actin binding, glutamate receptor activity, and cellular components such as focal adhesion, axon, and dendrite were enriched in both models. Nevertheless, more gene functions and pathways were significantly enriched in either of the two CSDS models but not the other. For example, biological processes related to cholesterol metabolism, neuron projection development, Golgi to vacuole transport and myelination were enriched in PS. In contrast, biological processes related to axonogenesis, and protein localization to membrane, and pathways related to GABAergic synapse and Retrograde endocannabinoid signaling were specifically enriched in ES.

FIGURE 9

Figure 9. Gene ontology and KEGG pathways associated with divergent G-to-A RNA editing patterns between PS and ES. The items with the most significant P-values are shown for (A) biological processes, (B) molecular functions, and (C) cellular components, as well as (D) KEGG pathways. ES, emotional stress; PS: physical stress.

Enrichment of Differential G-to-A RNA Editing in Genes With Multiple Editing Events

It was noted that a large proportion of genes in the mouse VTA RNA-Seq were observed to contain multiple G-to-A RNA editing events. The top ten genes with the largest counts of editing events are shown in Table 1. The microtubule associated protein 1b (Map1b) and HECT, UBA And WWE domain containing E3 ubiquitin protein ligase 1 (Huwe1) genes, were found to be the top two genes with the largest counts of editing events. Moreover, differential G-to-A RNA editing was also observed in the two genes. As illustrated in Table 2, Fisher’s exact test implicated that differential G-to-A RNA editing between control and CSDS or between PS and ES were dramatically enriched in genes with more than one editing event (OR = 2.6, 95% CI = 1.5–4.5, P = 2.26 × 10^–4 and OR = 2.6, 95% CI = 1.6–4.2, P = 5.08 × 10^–5).

TABLE 1

Table 1. Top ten genes with the largest counts of G-to-A RNA editing events.

TABLE 2

Table 2. Enrichment of differential G-to-A RNA editing in genes with multiple editing events.

Discussion

The rodent model of CSDS has been widely used in the study of stress-related diseases such as major depression and PTSD, providing an important and useful tool for the understanding of these diseases. The current study for the first time systematically investigated G-to-A RNA editing in the brain VTA in CSDS mice at the transcriptome-level, and demonstrated dynamic changes of G-to-A RNA editing associated with CSDS as well as divergent editing patterns found under ES or PS conditions.

Ventral tegmental area is a key region of the dopaminergic brain reward system, and its role in mood and anxiety disorders has been implicated by emerging evidence (Baik, 2020). Therefore, our results demonstrate that widespread G-to-A RNA editing is found in the mouse VTA transcriptome, indicating its potentially important biological functions in the brain. Nevertheless, the actual biological functions of G-to-A RNA editing, and how it is related to neurological and mental disorders remain largely unclear. By far a number of targets have been reported for G-to-A RNA editing, such as HNRNPK in colorectal cancer (Klimek-Tomczak et al., 2006) and WT1 in cord blood samples (Niavarani et al., 2015). For neural tissue and central nerve system, G-to-A editing was reported in GRIK3 in the fetal human brain (Nutt et al., 1994). Moreover, G-to-A editing led to a missense variant p.R441H in TPH2 (Grohmann et al., 2010) with decreased its enzyme activity, which had been associated with major depression disorder (Zhang et al., 2005). In the current study, multiple G-to-A editing events in the 3′-UTR of mouse Grik3 were observed. In addition, a high proportion of missense and stop-gain events were detected in genes with high expression levels in VTA. It is also worth noting that coding variants accounted for a high proportion (59.9%) of the RNA editing variants found, while missense, stop-gain and stop-lost variants accounted for another 68.2% of the coding variants. This may imply that G-to-A RNA editing could have a potential impact on the sequences of protein expressed in VTA. Such findings thus underscore the potential role of G-to-A editing in the brain.

Recently, RNA editing has been reported to be involved in CSDS. For example, by using target sequencing, Dick et al. (2019) reported dynamic regulation of A-to-I editing in 5-Hydroxytryptamine Receptor 2C (Htr2c) in CSDS mice. In line with such existing evidence, our results revealed changes of G-to-A editing events in CSDS mice, suggesting that CSDS could have dramatic effects on G-to-A editing in the brain. Moreover, divergent G-to-A editing patterns contributed to the difference between CSDS paradigms induced under ES and PS conditions. As a result, gene list enrichment further suggested that such changes may lead to altered gene functions or pathways in the brain.

It was noteworthy that a number of the genes with differential G-to-A editing in CSDS had been reported to be associated with neurological or psychiatric diseases, or functions as key components in neurons (Supplementary Table 2). In the current study, the most significant difference between CSDS and controls was observed in Fat1. It encodes a member of the cadherin superfamily. Two missense events in Fat1 were exclusive found in controls but not CSDS (Figure 4A and Supplementary Table 2), which might affect its protein function or structure., and its human ortholog FAT1 is associated with nicotine dependence and major depression (Blair et al., 2006). Likewise, the gene with the most significantly upregulated G-to-A RNA editing was found of Slc6a1, which encodes for the GABA transporter protein type 1 (also known as Gat1). The Slc6a1 transporter removes GABA, the primary inhibitory neurotransmitter in the CNS, from the synaptic cleft, and maintains low extracellular levels of GABA. Mutations of SLC6A1 were reported in a variety of neurological and psychiatric disorders, myoclonic atonic epilepsy, intellectual disability (Johannesen et al., 2018, 1), autism (Wang et al., 2020), as well as schizophrenia (Rees et al., 2020). Association of Slc6a1 with mood disorders had also been implicated. Slc6a1 deficient mice demonstrated reduced depression and anxiety-like behaviors (Liu et al., 2007; Gong et al., 2015). Moreover, genetic polymorphisms of SLC6A1 were reported to be associated with anxiety (Thoeringer et al., 2009). In line with such existing evidence, a trend of increased expression of Slc6a1 was found in CSDS mice with marginal significance (P = 0.087) together with increased G-to-A RNA editing in the in the 3′-UTR. Such findings, taken together, pointed to a potential effect of G-to-A RNA editing on Slc6a1 expression and GABAergic system in CSDS.

In particular, our findings re-emphasize the potential importance of RNA editing glutamate receptors, which are linked to fast excitatory neurotransmission in the CNS. Ionotropic glutamate receptors consist of three groups according to their specific agonists, includingα-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors, N-methyl-D-aspartate (NMDA) receptors, and kainate receptors. By far, Non-NMDA glutamate receptors have been known to be important targets for A-to-I RNA editing. Interestingly, our results implicate that NMDA including glutamate ionotropic receptor NMDA type subunit 2b (Grin2b) and glutamate ionotropic receptor NMDA type subunit associated protein 1 (Grina) could be the two most significant glutamate ionotropic receptors subjected to G-to-A RNA editing in CSDS. In addition, G-to-A RNA editing was also found in kainate receptors and AMPA receptors (Supplementary Table 1). GO analysis showed that the function of glutamate receptor-related gene functions was enriched in genes with G-to-A RNA editing in CSDS, suggesting that glutamate receptor-related genes may also be important targets of G-to-A RNA editing, which could be related to stress response.

In addition, divergent RNA-editing patterns were observed between ES and PS in the current study (Supplementary Table 3). The most significantly differential editing was found in Cul4a. It was noteworthy that a SNP in CUL4A was found to be associated with bipolar disorder in a genome-wide association study (GWAS) (Stahl et al., 2019). Moreover, a recent analysis revealed CUL4A as one of the 112 hub genes in the complex network of GWAS-identified genes in bipolar disorder (Xie et al., 2017). Therefore, the differential G-to-A RNA editing in Cul4a were probably in line with a potential role of the gene in emotional stress and mood disorders.

G-to-A RNA editing was not evenly distributed among expressed genes in the transcriptome. A number of genes were found to contain more than one editing event, suggesting that these genes might be editing hotspots. Moreover, genes with multiple editing events showed a significantly higher possibility to be differentially edited in CSDS compared to controls. These genes are mostly involved in development and function of the nervous system, or associated with neurological and mental disorders (Supplementary Tables 2, 3). For example, the Map1b gene encodes a protein member of the microtubule-associated protein family involved in essential microtubule assembly in neurogenesis. Loss-of-function mutations of MAP1B were reported in patients with intellectual disability (Walters et al., 2018). Another example was Huwe1, which encodes an E3 ubiquitin ligase, and plays a key role in cerebral cortex neurogenesis (Zhao et al., 2009; D’Arca et al., 2010). Mutations in this gene were reported in including both non-syndromic and syndromic forms of X-linked intellectual disability (Moortgat et al., 2018; Muthusamy et al., 2020).

In conclusion, the current study identified high-confidence G-to-A RNA editing in the VTA transcriptome of adult male mouse brains, and revealed dynamic changes in different CSDS models. With a number of stress/mood disorder-related genes revealed to be involved in altered G-to-A RNA editing, especially those with multiple editing events, our findings thus provide new insight into the understanding of the role played by brain G-to-A RNA editing in CSDS and its related diseases.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics Statement

The animal study was reviewed and approved by the Animal Ethics Committee of Jiangnan University.

Author Contributions

JT and C-YR performed the bioinformatic analysis. Z-YW improved the data analysis pipeline. JX and J-HC conceived the project and planned the experiments. FZ participated in the discussion. All authors contributed to the final manuscript.

Funding

This study was supported in part by grants from the National Natural Science Foundation of China (No. 31671311), the National first-class discipline program of Light Industry Technology and Engineering (LITE2018-14), the “Six Talent Peak” Plan of Jiangsu Province (No. SWYY-127), the Innovative and Entrepreneurial Talents of Jiangsu Province, the Program for High-Level Entrepreneurial and Innovative Talents of Jiangsu Province, Natural Science Foundation of Guangdong Province/Guangdong Basic and Applied Basic Research Foundation (2019A1515012062), Taihu Lake Talent Plan, and Fundamental Research Funds for the Central Universities (JUSRP51712B and JUSRP1901XNC), Wuxi Institute of Translational Medicine, and Postgraduate Research and Practice Innovation Program of Jiangsu Province (KYCX20_1946).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.680548/full#supplementary-material

Footnotes

1. ^ https://www.ncbi.nlm.nih.gov/geo/
2. ^ https://www.ensembl.org/vep
3. ^ http://srv00.recas.ba.infn.it/atlas/index.html
4. ^ https://maayanlab.cloud/Enrichr/

References

Baik, J.-H. (2020). Stress and the dopaminergic reward system. Exp. Mol. Med. 52, 1879–1890. doi: 10.1038/s12276-020-00532-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Blair, I. P., Chetcuti, A. F., Badenhop, R. F., Scimone, A., Moses, M. J., Adams, L. J., et al. (2006). Positional cloning, association analysis and expression studies provide convergent evidence that the cadherin gene FAT contains a bipolar disorder susceptibility allele. Mol. Psychiatry 11, 372–383. doi: 10.1038/sj.mp.4001784

PubMed Abstract | CrossRef Full Text | Google Scholar

Breen, M. S., Dobbyn, A., Li, Q., Roussos, P., Hoffman, G. E., Stahl, E., et al. (2019). Global landscape and genetic regulation of RNA editing in cortical samples from individuals with schizophrenia. Nat. Neurosci. 22, 1402–1412. doi: 10.1038/s41593-019-0463-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Christofi, T., and Zaravinos, A. (2019). RNA editing in the forefront of epitranscriptomics and human health. J. Transl. Med. 17:319. doi: 10.1186/s12967-019-2071-4

PubMed Abstract | CrossRef Full Text | Google Scholar

D’Arca, D., Zhao, X., Xu, W., Ramirez-Martinez, N. C., Iavarone, A., and Lasorella, A. (2010). Huwe1 ubiquitin ligase is essential to synchronize neuronal and glial differentiation in the developing cerebellum. Proc. Natl. Acad. Sci. U.S.A. 107, 5875–5880. doi: 10.1073/pnas.0912874107

PubMed Abstract | CrossRef Full Text | Google Scholar

Dick, A. L. W., Khermesh, K., Paul, E., Stamp, F., Levanon, E. Y., and Chen, A. (2019). Adenosine-to-Inosine RNA editing within corticolimbic brain regions is regulated in response to chronic social defeat stress in mice. Front. Psychiatry 10:277. doi: 10.3389/fpsyt.2019.00277

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., and Gingeras, T. R. (2015). Mapping RNA-seq reads with STAR. Curr. Protoc. Bioinform. 51, 11.14.1–11.14.19. doi: 10.1002/0471250953.bi1114s51

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, X., Shao, Y., Li, B., Chen, L., Wang, C., and Chen, Y. (2015). γ-aminobutyric acid transporter-1 is involved in anxiety-like behaviors and cognitive function in knockout mice. Exp. Ther. Med. 10, 653–658. doi: 10.3892/etm.2015.2577

PubMed Abstract | CrossRef Full Text | Google Scholar

Grohmann, M., Hammer, P., Walther, M., Paulmann, N., Büttner, A., Eisenmenger, W., et al. (2010). Alternative splicing and extensive rna editing of human TPH2 transcripts. PLoS One 5:e8956. doi: 10.1371/journal.pone.0008956

PubMed Abstract | CrossRef Full Text | Google Scholar

Gu, T., Buaas, F. W., Simons, A. K., Ackert-Bicknell, C. L., Braun, R. E., and Hibbs, M. A. (2012). Canonical A-to-I and C-to-U RNA editing is enriched at 3’UTRs and microRNA target sites in multiple mouse tissues. PLoS One 7:e33720. doi: 10.1371/journal.pone.0033720

PubMed Abstract | CrossRef Full Text | Google Scholar

Johannesen, K. M., Gardella, E., Linnankivi, T., Courage, C., Martin, A., de, S., et al. (2018). Defining the phenotypic spectrum of SLC6A1 mutations. Epilepsia 59, 389–402. doi: 10.1111/epi.13986

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawahara, Y., Ito, K., Sun, H., Aizawa, H., Kanazawa, I., and Kwak, S. (2004). RNA editing and death of motor neurons. Nature 427, 801–801. doi: 10.1038/427801a

PubMed Abstract | CrossRef Full Text | Google Scholar

Klimek-Tomczak, K., Mikula, M., Dzwonek, A., Paziewska, A., Karczmarski, J., Hennig, E., et al. (2006). Editing of hnRNP K protein mRNA in colorectal adenocarcinoma and surrounding mucosa. Br. J. Cancer 94, 586–592. doi: 10.1038/sj.bjc.6602938

PubMed Abstract | CrossRef Full Text | Google Scholar

Koboldt, D. C., Zhang, Q., Larson, D. E., Shen, D., McLellan, M. D., Lin, L., et al. (2012). VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576. doi: 10.1101/gr.129684.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuleshov, M. V., Jones, M. R., Rouillard, A. D., Fernandez, N. F., Duan, Q., Wang, Z., et al. (2016). Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 44, W90–W97. doi: 10.1093/nar/gkw377

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, G.-X., Cai, G.-Q., Cai, Y.-Q., Sheng, Z.-J., Jiang, J., Mei, Z., et al. (2007). Reduced anxiety and depression-like behaviors in mice lacking GABA transporter subtype 1. Neuropsychopharmacology 32, 1531–1539. doi: 10.1038/sj.npp.1301281

PubMed Abstract | CrossRef Full Text | Google Scholar

Mansi, L., Tangaro, M. A., Lo Giudice, C., Flati, T., Kopel, E., Schaffer, A. A., et al. (2021). REDIportal: millions of novel A-to-I RNA editing events from thousands of RNAseq experiments. Nucleic Acids Res. 49, D1012–D1019. doi: 10.1093/nar/gkaa916

PubMed Abstract | CrossRef Full Text | Google Scholar

McCarthy, D. J., Chen, Y., and Smyth, G. K. (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40, 4288–4297. doi: 10.1093/nar/gks042

PubMed Abstract | CrossRef Full Text | Google Scholar

McLaren, W., Gil, L., Hunt, S. E., Riat, H. S., Ritchie, G. R. S., Thormann, A., et al. (2016). The ensembl variant effect predictor. Genome Biol. 17:122. doi: 10.1186/s13059-016-0974-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Moortgat, S., Berland, S., Aukrust, I., Maystadt, I., Baker, L., Benoit, V., et al. (2018). HUWE1 variants cause dominant X-linked intellectual disability: a clinical study of 21 patients. Eur. J. Hum. Genet. 26, 64–74. doi: 10.1038/s41431-017-0038-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Muthusamy, B., Nguyen, T. T., Bandari, A. K., Basheer, S., Selvan, L. D. N., Chandel, D., et al. (2020). Exome sequencing reveals a novel splice site variant in HUWE1 gene in patients with suspected Say-Meyer syndrome. Eur. J. Med. Genet. 63:103635. doi: 10.1016/j.ejmg.2019.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Niavarani, A., Currie, E., Reyal, Y., Anjos-Afonso, F., Horswell, S., Griessinger, E., et al. (2015). APOBEC3A is implicated in a novel class of G-to-A mRNA editing in WT1 transcripts. PLoS One 10:e0120089. doi: 10.1371/journal.pone.0120089

PubMed Abstract | CrossRef Full Text | Google Scholar

Nutt, S. L., Hoo, K. H., Rampersad, V., Deverill, R. M., Elliott, C. E., Fletcher, E. J., et al. (1994). Molecular characterization of the human EAA5 (GluR7) receptor: a high-affinity kainate receptor with novel potential RNA editing sites. Recept. Channels 2, 315–326.

Google Scholar

Peña, C. J., Kronman, H. G., Walker, D. M., Cates, H. M., Bagot, R. C., Purushothaman, I., et al. (2017). Early life stress confers lifelong stress susceptibility in mice via ventral tegmental area OTX2. Science 356, 1185–1188. doi: 10.1126/science.aan4491

PubMed Abstract | CrossRef Full Text | Google Scholar

Peña, C. J., Smith, M., Ramakrishnan, A., Cates, H. M., Bagot, R. C., Kronman, H. G., et al. (2019). Early life stress alters transcriptomic patterning across reward circuitry in male and female mice. Nat. Commun. 10:5098. doi: 10.1038/s41467-019-13085-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Rees, E., Han, J., Morgan, J., Carrera, N., Escott-Price, V., Pocklington, A. J., et al. (2020). De novo mutations identified by exome sequencing implicate rare missense variants in SLC6A1 in schizophrenia. Nat. Neurosci. 23, 179–184. doi: 10.1038/s41593-019-0565-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Small, I. D., Schallenberg-Rüdinger, M., Takenaka, M., Mireau, H., and Ostersetzer-Biran, O. (2020). Plant organellar RNA editing: what 30 years of research has revealed. Plant J. 101, 1040–1056. doi: 10.1111/tpj.14578

PubMed Abstract | CrossRef Full Text | Google Scholar

Stahl, E. A., Breen, G., Forstner, A. J., McQuillin, A., Ripke, S., Trubetskoy, V., et al. (2019). Genome-wide association study identifies 30 loci associated with bipolar disorder. Nat. Genet. 51, 793–803.

Google Scholar

Tan, M. H., Li, Q., Shanmugam, R., Piskol, R., Kohler, J., Young, A. N., et al. (2017). Dynamic landscape and regulation of RNA editing in mammals. Nature 550, 249–254. doi: 10.1038/nature24041

PubMed Abstract | CrossRef Full Text | Google Scholar

Thoeringer, C. K., Ripke, S., Unschuld, P. G., Lucae, S., Ising, M., Bettecken, T., et al. (2009). The GABA transporter 1 (SLC6A1): a novel candidate gene for anxiety disorders. J. Neural. Transm. 116, 649–657. doi: 10.1007/s00702-008-0075-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Van der Auwera, G. A., Carneiro, M. O., Hartl, C., Poplin, R., del Angel, G., Levy-Moonshine, A., et al. (2013). From FastQ data to high confidence variant calls: the genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinform. 11, 11.10.1–11.10.33. doi: 10.1002/0471250953.bi1110s43

PubMed Abstract | CrossRef Full Text | Google Scholar

Walters, G. B., Gustafsson, O., Sveinbjornsson, G., Eiriksdottir, V. K., Agustsdottir, A. B., Jonsdottir, G. A., et al. (2018). MAP1B mutations cause intellectual disability and extensive white matter deficit. Nat. Commun. 9:3456. doi: 10.1038/s41467-018-05595-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Poliquin, S., Mermer, F., Eissman, J., Delpire, E., Wang, J., et al. (2020). Endoplasmic reticulum retention and degradation of a mutation in SLC6A1 associated with epilepsy and autism. Mol. Brain 13:76.

Google Scholar

Warren, B. L., Vialou, V. F., Iñiguez, S. D., Alcantara, L. F., Wright, K. N., Feng, J., et al. (2013). Neurobiological sequelae of witnessing stressful events in adult mice. Biol. Psychiatry 73, 7–14. doi: 10.1016/j.biopsych.2012.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Z., Yang, X., Deng, X., Ma, M., and Shu, K. (2017). A genome-wide association study and complex network identify four core hub genes in bipolar disorder. Int. J. Mol. Sci. 18:2763.

Google Scholar

Zhang, X., Gainetdinov, R. R., Beaulieu, J.-M., Sotnikova, T. D., Burch, L. H., Williams, R. B., et al. (2005). Loss-of-function mutation in tryptophan hydroxylase-2 identified in unipolar major depression. Neuron 45, 11–16. doi: 10.1016/j.neuron.2004.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Arca, D. D., Lim, W. K., Brahmachary, M., Carro, M. S., Ludwig, T., et al. (2009). The N-Myc-DLL3 cascade is suppressed by the ubiquitin ligase Huwe1 to inhibit proliferation and promote neurogenesis in the developing brain. Dev. Cell 17, 210–221. doi: 10.1016/j.devcel.2009.07.009

PubMed Abstract | CrossRef Full Text | Google Scholar