Adenosine-to-inosine (A-to-I) editing is hypothesized to facilitate adaptive evolution by expanding proteomic diversity through an epigenetic approach. However, it is challenging to provide evidences to support this hypothesis at the whole editome level. In this study, we systematically characterized 2,114 A-to-I RNA editing sites in female and male brains of D. melanogaster, and nearly half of these sites had events evolutionarily conserved across Drosophila species. We detected strong signatures of positive selection on the nonsynonymous editing sites in Drosophila brains, and the beneficial editing sites were significantly enriched in genes related to chemical and electrical neurotransmission. The signal of adaptation was even more pronounced for the editing sites located in X chromosome or for those commonly observed across Drosophila species. We identified a set of gene candidates (termed “PSEB” genes) that had nonsynonymous editing events favored by natural selection. We presented evidence that editing preferentially increased mutation sequence space of evolutionarily conserved genes, which supported the adaptive evolution hypothesis of editing. We found prevalent nonsynonymous editing sites that were favored by natural selection in female and male adults from five strains of D. melanogaster. We showed that temperature played a more important role than gender effect in shaping the editing levels, although the effect of temperature is relatively weaker compared to that of species effect. We also explored the relevant factors that shape the selective patterns of the global editomes. Altogether we demonstrated that abundant nonsynonymous editing sites in Drosophila brains were adaptive and maintained by natural selection during evolution. Our results shed new light on the evolutionary principles and functional consequences of RNA editing.
Genomic mutations are the major sources for phenotypic changes and adaptation [1–4]. In diploid multicellular organisms, a nonsynonymous DNA mutation (a mutation that alters the amino acid sequence of a protein) will permanently affect the protein products in all the cells (soma or germline) that express the mutant allele. The “all-or-none” property of DNA mutations might incur pleiotropic effects that are antagonistic among cell types, tissues, developmental stages, sexes, or other aspects of life history [5–7], which would constrain the available genetic diversity for a species. Given the prevalence of pleiotropic effects in the genomes [8–10], the sequence space might be inaccessible to many mutations, which potentially slows down the rate of phenotypic evolution and adaptation [11]. However, the transcriptomic or proteomic diversity limited by mutation sequence space could be expanded by the alteration of RNA sequences in an epigenetic approach, such as RNA editing, which was hypothesized to facilitate adaptation [12–14]. In addition, RNA editing has the advantage to quickly respond to environmental stress and adjust the activity of final protein products accordingly [15,16]. RNA editing is an evolutionarily conserved mechanism that alters RNA sequences at the co-transcriptional or post-transcriptional level [13,17–23]. Among various RNA editing systems in animals, the base substitution from adenosine (A) to inosine (I), termed A-to-I editing, is the most common form [13,20]. Due to the high level of structural similarity between inosine (I) and guanosine (G), the cellular machineries, such as ribosomes, spliceosomes or the microRNA ribonucleoprotein complex (miRNP), would recognize I as G during translation [13,20,24,25], splicing [26–29], microRNA target recognition [30–32], or other RNA biological processes [14]. Therefore, A-to-I RNA editing usually produces a change similar to an A-to-G DNA change in particular tissues or developmental stages, which potentially increases phenotypic plasticity without the alteration of genomic sequences [13,20,24,25]. The adenosine deaminase acting on RNA (ADAR) family are the enzymes that convert adenosine (A) to inosine (I) in pre-mRNAs [33–36]. Although multiple Adar genes are encoded in the genomes of mammals and worms, there is only one Adar locus in Drosophila [37,38], which is predominately expressed in the nervous system [39]. The substrates of ADAR are usually double-stranded RNAs [34,36,40–42]. A-to-I editing plays essential roles in many biological processes [18,19,43–45], and the abolition of Adar in D. melanogaster severely affects its viability and behavior [33,34,46]. Previous studies have identified thousands of A-to-I editing sites in different developmental stages, adult heads or whole animals of Drosophila [47–52]. In addition, A-to-I editing has been extensively characterized in other organisms, such as humans [53–58], macaques [59,60], mice [61], worms [62], and squids [63]. Despite these intriguing advances, only a few examples of the advantageous effects conferred by RNA editing have been demonstrated [13,14,20,63]. For example, the A-to-I editing events in Kv1 mRNA provide numerous adaptive amino acid changes that allow the octopus to adapt to extremely cold temperatures [64]. The functional consequences of the majority of A-to-I editing events, however, remain to be explored. In fact, comparative genomics has demonstrated that only a small fraction of the human A-to-I editing events were evolutionarily conserved [65–68]. Furthermore, it was nicely demonstrated that the editing events in primate coding regions were generally non-adaptive [60,67,68]. Nevertheless, the targets of RNA editing might have evolved rapidly across species because A-to-I editing in mammals predominantly occurs in repetitive sequences [53–55,61,65,66], while the editing events in Drosophila are mainly located in coding regions of genes encoding neurotransmitters or ion channels [47–50,69,70]. Therefore, the evolutionary forces acting on A-to-I RNA editing might be different between Drosophila [50] and primates [60,67,68]. If A-to-I editing indeed facilitates adaptation by expanding proteomic diversity, we expect to observe predominant signals of adaptation in the editing sites. A recent study [50] reported that although signals of positive selection could be found in genes of the nervous system, the A-to-I RNA editing events were overall subject to purifying selection in Drosophila. Additionally, the overall effect of natural selection on the editome is different across Drosophila developmental stages [50]. Despite these intriguing discoveries, it still remains a mystery whether or not we can find evidence to show that the whole editome is overall adaptive. Specifically, we are interested in testing whether the nonsynonymous A-to-I editing events in Drosophila brains, the core component of the nervous system, are predominantly adaptive. Furthermore, several other fundamental questions on editing deserve to be further investigated: 1) Do editing sites preferentially increase sequence space of evolutionarily conserved genes? 2) Why does the global editome of different tissues or developmental stages show differential selective patterns? 3) How does temperature shape the global editomes? Answers to these questions will help understand the evolutionary principles and functional consequences of RNA editing. In this study we addressed these questions by systematically sequencing the transcriptomes and deciphering A-to-I editing in the female and male brains of three Drosophila species at different temperatures. With evolutionary analysis from different perspectives, we provided lines of evidence to demonstrate that the nonsynonymous editing events in coding regions are generally adaptive in brains of Drosophila. Then we identified a set of gene candidates that had nonsynonymous editing events favored by natural selection. Overall our results demonstrated that abundant nonsynonymous editing events in Drosophila brains were adaptive and maintained by natural selection during evolution.
To comprehensively characterize the A-to-I editing landscapes in brains of Drosophila, we set out to dissect the brains of 1- to 5-day-old or 1- to 14-day-old female and male adults of the inbred ISO-1 strain of Drosophila melanogaster that were constantly raised at 25°C, or raised at 25°C and treated at 30°C for 14 hours or 48 hours (Table 1). Next we selected the poly (A) -tailed mRNAs, fragmented them, ligated the mRNA fragments with adaptors, and deep sequenced the transcriptome of each brain sample (Materials and Methods). We obtained 13. 9–21. 6M reads mapped on the reference genome (see Table 1 and S1 Table for detailed statistics), and the median coverage on an exonic site in a library ranges from 5 to 9 reads (S1 Fig). As justified previously [71], the mRNA fragmentation library preparation procedure we employ minimizes the bias of non-uniform sequencing read coverage along mRNAs, which would reduce the bias in detecting editing in 3' ends of mRNAs. It is a challenging task to reliably distinguish the A-to-I editing events from SNPs [72–76], therefore, the ISO-1 strain used in this study, which was inbred and sequenced to assemble the reference genome of D. melanogaster [77], enables us to detect DNA-RNA differences with high accuracy and sensitivity. We employed a two-step strategy to identify editing sites in the brains of D. melanogaster (Fig 1A). First, we used the GATK RNA-Seq variant calling pipeline [78] to identify the candidate A-to-I editing sites in each brain library (i. e. , the A-to-G differences in the final sequencing results). We identified 1,531 unique sites with A-to-G DNA-RNA differences in these brain libraries, and such differences accounted for 81. 5% (with a standard error of 0. 97%) of all the differences detected by the GATK pipeline in each library (S2A Fig). In contrast, the proportion of A-to-G DNA differences (reference vs. alternative allele) was only 9. 9% out of all the mutations (S2B Fig) in the global populations of D. melanogaster [79]. This comparison justified the reliability and accuracy of our procedures in defining the candidate A-to-I editing sites. Second, we retrieved a total of 5,389 editing sites characterized in D. melanogaster in previous studies (972 in Graveley et al. [47], 1,350 in Rodriguez et al. [49], 3,581 in St. Laurant et al. [48], and 1,298 in Yu et al. [50]). Altogether, we obtained 5,925 unique candidate sites (986 sites overlapped between GATK and the four previous studies, S2 Table). For each candidate site in a brain library k, we calculate the probability that the A-to-G difference (if detected) is caused by editing with Pk (E1) = 1 − Pk (E0), where Pk (E0) is the probability that the difference is solely caused by sequencing error (ε), by incorporating the sequencing coverage (Ck) and the number of G allele (Lk) at that site. Next, we utilize the multiple-sample information and calculate the joint probability that this site is edited in at least one library, P (E1) = 1 − P (E0), where P (E0) is the probability that the A-to-G differences observed in that site across all the applicable libraries are entirely caused by sequencing errors (Materials and Methods). Among the 5,925 candidate sites, we did not detect expression of 664 sites in any brain library. For the remaining 5,261 expressed sites, we divided them into five exclusive classes with decreasing confidence based on P (E1), sequencing coverage, and the number of libraries in which the editing events were detected. Briefly, Class I (1,702 sites) were defined with the following criteria: 1) at FDR of 0. 001,2) the maximum sequencing coverage across all the libraries (Cmax) ≥ 10,3) the total coverage across all the libraries (Ctotal) ≥ 40, and 4) editing was detected in at least two libraries. Among the remaining sites that had editing events detected in at least two libraries, we defined Class II (447 editing sites) with these criteria: 1) at FDR of 0. 01,2) Cmax ≥ 5 and 3) Ctotal ≥ 16 (we also employed other cutoffs to define editing sites in Class I and II, and obtained results not very different from the results reported here; see S3 Table for details). 824 sites do not meet the aforementioned two criteria but have P (E1) > 0. 99, which suggests they might also be edited, although with lower confidence in brains of D. melanogaster (Class III). Moreover, 131 sites have editing detected in at least one library but have P (E1) ≤ 0. 99 (Class IV). Notably, we detected mRNA-Seq reads covering the remaining 2,157 candidate sites but none of them has editing events detected (Class V). The detailed information about these sites is presented in S2 Table. It is not surprising that the sequencing coverage decreases in the order of Class I, II, and III in each library [the median coverage in each library is 17. 4±1. 4 (mean ± s. e. throughout this study), 3. 87±0. 35 and 1±0 raw reads, respectively; S3A Fig]. Interestingly, although the sites in Class II have significantly lower coverage compared to sites in Class I (P < 0. 05 in each library, Kolmogorov-Smirnov tests), the editing levels are even significantly higher in Class II than in Class I (the median editing level in each library is 0. 24±0. 02 vs. 0. 16±0. 01 in Class II vs. I, S3B Fig). From another perspective, 45. 5% of the Class I sites were edited in all the eight brain libraries, meanwhile, only 8. 5% of the Class II sites were edited in all the eight brain libraries (P < 0. 01, Fisher’s exact test; S3C Fig). Compared to Class I and II, sites in Class III have both lower coverage and editing levels (S3A and S3B Fig). Sites in Class IV are extremely lowly edited and sites in Class V do not have any editing event detected in our samples; however, these two classes do not have the lowest sequencing coverage compared to the other three classes (the median coverage in a library is 18. 6±1. 1 and 6. 1±0. 35 for Class IV and V respectively, S3A Fig), suggesting they might have negligible editing in brains of D. melanogaster. To estimate the false positive rates of the editing sites in each class, we analyzed the RNA-Seq datasets from paired samples of wild-type strain w1118 vs. Adar5G1 mutant of D. melanogaster as conducted previously [50,51]. We found 1,145,161, and 103 editing sites in Class I, II, and III respectively that have editing events detected in w1118 heads, and correspondingly, 33,2, and 7 of these sites were detected in the heads of Adar5G1 mutant, yielding a false positive rate of 2. 88%, 1. 24% and 6. 80% for Class I, II, and III, respectively. Therefore, the sites in Class I and II captured the editing events in brains of D. melanogaster with high accuracy, and Class III sites were not considered in the down-stream analysis due to the high positive rate. For the sites in Class I and II, we identified 1,630 (1,243 in Class I and 387 in Class II) sites overlapped with previous studies [47–50], and 519 sites (459 in Class I and 60 in Class II) are novel in this study (S2 Table). It is not uncommon that many editing sites are not overlapping between studies in Drosophila [47–50]: on average 30. 4±3. 6% of the editing sites are shared in pairwise comparisons (ranging from 12. 8% to 54. 7%, S4 Table); and we observed comparable proportions of shared sites between our study and the previous ones: 21. 8%, 36. 7%, 61. 1% and 28. 2% of the Class I+II sites in our study are overlapping with Graveley et al. [47], Rodriguez et al. [49], St. Laurant et al. [48], and Yu et al. [50], respectively (S4 Table). Importantly, when we pooled Class I and II together, we found the novel sites have comparable false positive rates as the common ones in the w1118 vs. Adar5G1 mutant analysis (8/242 = 3. 31% vs. 27/1064 = 2. 54% for the novel vs. common sites). Furthermore, 111 of the novel sites are annotated in Ramaswami et al. [52], which is independent from this study. Taken together, we identified 2,114 “high-confidence” editing sites after combining sites in Class I and II (35 sites that have A-to-G difference in Adar5G1 mutants were removed), including 1,603 (75. 8%) sites overlapped with sites identified by previous studies [47–50] and 511 (24. 2%) novel sites (Fig 1B). The novel sites have slightly higher sequencing coverage than the common sites in the brain libraries (the median coverage is 25. 1±1. 5 and 20. 5±1. 1 for the former and latter, respectively, P < 0. 05 in each library, KS tests; S3D Fig), but generally lower editing levels (the median is 0. 18±0. 006 vs. 0. 36±0. 006, P < 10−16 in each library, KS tests; Fig 1C). Moreover, compared to the common editing sites, the novel sites are generally edited in fewer brain samples: 42. 9% of the common ones were detected in all the eight brain libraries, while only 20. 0% of the novel sites were detected in all the libraries (P < 10−16, Fisher’s exact test; S3E and S3F Fig). Altogether these results suggest that these novel editing sites are genuine but lowly edited in the brains, and were probably diluted in the samples of previous studies that were carried out in heads or whole flies [47–50]. Among the 2,114 high-confidence sites (Fig 1D), 235 (11. 1%) are in intergenic regions, 42 (2. 0%) are in ncRNAs, and 1,837 (86. 9%) are in 517 protein-coding genes, including 20 (0. 95%) in 5' UTRs, 550 (26. 0%) in introns, 414 (19. 6%) in 3' UTRs, 678 (32. 1%) nonsynonymous (in CDS regions and causes amino acid changes when edited, abbreviated as N throughout this study), and 144 (6. 8%) synonymous (in CDS regions but do not cause amino acid changes when edited, abbreviated as S), and one editing site (chr3R: 18806029) that putatively disrupts the stop codon of CG18208 (UAG>UGG). The detailed annotation in each library was presented in S5 Table. The gene ontology analysis revealed that the high-confidence exonic editing sites were significantly enriched in genes that encode transporters, synaptic vesicles or neurotransmitters (S6 Table and S7 Table for male and female brains, respectively; and the top 50 genes that had the largest number of editing sites were presented in S8 Table). For the exonic editing sites, the editing levels (averaged across libraries) decrease in the order of N (0. 319±0. 010), S (0. 214±0. 017), 3' UTRs (0. 168±0. 008), and 5' UTRs (0. 133±0. 020), with levels in N sites significantly higher than those in the other three categories in the brains of D. melanogaster (P < 0. 001, KS test; Fig 1D), suggesting that high levels of nonsynonymous editing events are overall favorable. Among the 550 intronic editing sites, 167 of them might also be exonic due to alternative splicing (we only used annotations of the canonical transcript and some intronic sites in the canonical transcripts might be coding in the non-canonical transcripts), and the coverage and editing levels (0. 336±0. 013) are comparable to the N sites (Fig 1D). Interestingly, the remaining 383 authentic intronic sites generally have significantly lower coverage than the coding regions (Fig 1D), however, high editing levels in these sites (0. 418±0. 006) were observed, supporting previous results that editing is exerted co-transcriptionally [49]. We uncovered a tendency that A-to-I editing events were more readily detected in the genes with higher expression levels (or adenosine sites with higher mRNA-Seq coverage). In each brain sample, when we grouped the expressed genes into 20 bins with increasing expression levels (only genes with RPKM ≥ 1 were considered), we found a significant positive correlation between the editing density (hereafter defined as the number of edited out of the total adenosine sites) and the gene expression level (P < 0. 001 in each library; S4A Fig). Similar patterns were observed if we grouped all the adenosine sites with increasing mRNA-Seq coverage in each sample (only sites ≥ 5X coverage were considered; S4B Fig). Analogous results were obtained if we weighted each editing site with its editing level (“level-weighted density of editing sites”, see S4A and S4B Fig). Consistent with previous observations [47–49], we found the editing density was significantly increased from the 5' to 3' of pre-mRNAs. After dividing the adenosine sites (≥ 5X coverage) into 20 equal bins along their positions in pre-mRNAs, our meta-gene analysis indicated that the editing density in each bin was significantly positively correlated with the relative distance of that bin to the transcriptional start sites (P < 0. 005 in each library; S4C Fig). Despite our experimental optimization, the poly (A) selection procedure still caused slightly increased coverage bias towards 3' ends of mRNAs (S4D Fig). However, we found the coverage difference between 5' and 3' of mRNAs was not the main cause of elevated editing density in the 3' ends of mRNAs with two analyses. First, in each library, we split each gene (RPKM ≥ 1) into two equal parts, calculated the RPKM values for each half-gene separately, ranked all the half-genes with increasing RPKM values, and grouped them into 20 bins. Next, in each bin, we combined the 5' and 3' half-genes independently and calculated the editing density in the 5' and 3' half. We found within each bin the editing density in the 3' half-genes are significantly higher than the 5' half genes (P < 0. 05 in each library; paired t tests, S4E Fig). To further reduce the coverage variation within the half-genes, we ranked all the adenosine sites (≥ 5X) with increasing coverage and binned them into 20 groups, and in each group we calculated the editing density in the 5' (front) half and 3' (rear) half of pre-mRNAs independently. We also found the editing density were significantly higher for sites in the rear half compared to sites in the front half of pre-mRNAs (P < 0. 001 in each library; paired t tests, S4F Fig). Taken together, the increasing editing density along mRNAs is not likely caused by detection bias, but more likely shaped by the recruitment of ADAR to the transcription elongation complex, as previous functional studies demonstrated [49,80]. We predicted that 591 (50. 1%) of the 1,179 exonic editing sites were located in stable local hairpin structures of mRNAs (Materials and Methods), such as Adar (S5A Fig), rtp, DIP1, rdgA, CG43897, and CG42540 (editing events in these genes were verified with Sanger sequencing; S5B and S5C Fig). In contrast, we obtained only 363 exonic sites (332–393 sites within 95% CI) located in stable hairpin structures after comprehensively folding all the transcripts expressed in brains and randomly sampling the equal amounts of editing sites (Materials and Methods; Fig 1E). Similar results were obtained when we focused on the N or S editing sites individually (P < 0. 002 in simulations for both cases; Fig 1E). In addition, we found 181 intronic editing sites located in stable hairpin structures when we folded the pre-mRNA sequences. Long-range pseudoknots are another class of RNA substrates recognized by ADAR [81]. By extensively folding the flanking sequences of the editing sites (see Methods for details), we inferred 260 (22. 1%) exonic editing sites that were located outside stable hairpin structures but were located in stems of long-range pseudoknots in pre-mRNAs of genes, such as the 3' UTR of Adar (S5D Fig), nrm, B52, nAchRbeta1, CG8034 and roX1 (S6 Fig; the editing events in nrm were verified by Sanger sequencing of the cDNA and genomic DNA, S6B and S6C Fig). Taken together, our results systematically demonstrated that at least 874 (74. 1%) of the exonic A-to-I editing sites in the brains of D. melanogaster were located in pre-mRNA regions that formed stable secondary structures. These results also well explain why the A-to-I editing sites are located in clusters, as commonly observed in previous studies [41,47–49,62]. By clustering the editing sites with distances smaller than 100 nucleotides, we identified a total of 1,320 editing sites that form 413 clusters in brains of D. melanogaster (S7 Fig), and unusually large editing clusters were frequently observed, such as in NaCP60E and CaMKII (the Sanger verification was presented in S8 Fig). To characterize the A-to-I editing events that were evolutionarily conserved (i. e. , commonly observed) across species, we deep sequenced the poly (A) -tailed transcriptomes of female and male brains of 1- to 5-day-old D. simulans and D. pseudoobscura that were accommodated at the same temperature conditions as D. melanogaster (six libraries for each species). The mapped reads range from 8. 7–16. 4M in each library of D. simulans, and 10. 9–17. 8M in D. pseudoobscura (Table 2, Table 3 and S1 Table for detailed information), and the median sequencing coverage on an exonic site in a library ranges from 5 to 9 reads in D. simulans, and ranges from 4 to 5 in D. pseudoobscura (S1 Fig). D. simulans diverged from D. melanogaster ~5. 4 million years ago (Fig 2) while D. pseudoobscura diverged from D. melanogaster approximately 55 million years ago [82]. Comparing A-to-I editing across these three species will help us understand the role of natural selection in shaping the brain editomes during evolution. To exclude SNPs in the RNA editing characterization, we also deep sequenced the genomic DNA of the same strain of D. simulans (the median coverage per site is 46, totally 313,133 SNPs, S9A Fig) and D. pseudoobscura (the median coverage per site is 47, totally 489,828 SNPs, S9B Fig) and masked all the SNPs (Materials and Methods). For each high-confidence editing site in brains of D. melanogaster, we employed two complementary approaches to search for the orthologous sites in D. simulans and D. pseudoobscura. First, we used liftOver [83] to convert the genomic coordinates of the orthologous sites between D. melanogaster and D. simulans, or between D. melanogaster and D. pseudoobscura, based on the pairwise genome alignments as previously conducted [51] (termed “g_align” approach, Materials and Methods). Second, we parsed out the genomic coordinates with the pairwise CDS alignments that were made based on the protein alignments between D. melanogaster and the other species (“c_align” approach). We pooled orthologous sites by the two approaches together. For each site in each species, we calculated the joint probability that this site is edited in at least one library P (E1). At FDR of 0. 05, we identified 996 sites edited in D. simulans (S9 Table), and 451 sites edited in D. pseudoobscura (S10 Table), and 367 sites edited in both D. simulans and D. pseudoobscura (Fig 2). We present the editing sites evolutionarily conserved in the same gender under the same temperature conditions in D. simulans (Table 2 and S11 Table) and D. pseudoobscura (Table 3 and S12 Table). For the editing sites we characterized in the brains of D. melanogaster, 34. 3–44. 0% of them have editing events detected in the matched samples of D. simulans (Table 2), and 22. 3–24. 1% of them have editing in the matched samples of D. pseudoobscura (Table 3). Note the proportion of editing sites in D. melanogaster that have editing events detected in brains of other species varies across libraries since we required the sites are edited in both paired samples. In general, with divergence increases, the level of conserved editing sites decreased, suggesting the editing events are evolutionary dynamic. Notably, for the 996 editing sites with conserved events in both D. melanogaster and D. simulans, and the 451 editing sites with conserved events in both D. melanogaster and D. pseudoobscura, we found 416 (41. 8%) and 78 (17. 3%) of them are located outside the coding regions, respectively (Tables 2 and 3), which is consistent with a recent study [84] and suggests a possible functional role for these sites, such as influencing alternative splicing [26–29], microRNA targeting [30–32], or other cellular processes related to RNAs [14]. Comparing the editing sites with conserved and non-conserved events revealed two interesting features. First, in each brain library, the editing levels are significantly higher in the sites with evolutionarily conserved editing events than in the remaining sites (the mean level in a library is 0. 340±0. 008 vs. 0. 252 ± 0. 009 in the D. melanogaster/D. simulans comparison, and 0. 323±0. 012 vs. 0. 187±0. 009 in the D. melanogaster/D. pseudoobscura comparison; P < 0. 01 in each comparison, KS tests, S10 Fig). Second, the N sites are significantly enriched in the editing sites that are evolutionarily conserved: 72. 9% (494 out of 678) N sites compared to 35. 0% (502 out of 1436) of the remaining sites that have evolutionarily conserved events between D. melanogaster and D. simulans (P < 10−10, Fisher’s exact test), and 47. 9% (325 out of 678) N sites compared to 8. 77% (126 out of 1436) of the remaining sites that have evolutionarily conserved events between D. melanogaster and D. pseudoobscura (P < 10−10, Fisher’s exact test), suggesting the nonsynonymous editing events are generally maintained and regulated by different evolutionary forces compared to the other sites. There are 84 editing sites that have editing events detected in both D. melanogaster and D. pseudoobscura but without editing events confidently identified in D. simulans. Nevertheless, this does not necessarily mean these sites are not edited in D. simulans (for 60 of these sites we did not find the orthologous sites in D. simulans, and for the 24 remaining sites, 10 of them have low level of editing but are undistinguishable from sequencing errors; S13 Table). Sampling bias frequently causes the sites with low expression or low editing levels to yield no editing signals in the sequencing libraries. Therefore, next we only focused on the sites with high sequencing coverage to explore the possible gain and loss patterns of editing events. We obtained 87 editing sites that have minimal editing level of 0. 05 in D. melanogaster and have at least 200 raw reads (across all the libraries) in both D. simulans and D. pseudoobscura. We found 52 sites with editing events reliably detected in all the three species. For each of the remaining 35 sites, in case no editing event was discovered at a site in a sample m in D. simulans (or D. pseudoobscura), we calculate Pm (D0), the probability that this observation happens by sampling bias or because the editing signal was abolished by sequencing error (ε), given a depth of Cm and an assumed editing level lm at that site (Materials and Methods). We assumed the orthologous sites in the other species have the same editing level as in D. melanogaster and calculated the joint probability P (D0) that a site was edited despite zero edited allele was detected in all the libraries. After correcting for multiple testing, at FDR of 0. 05, we found 20 sites with editing present in both D. melanogaster and D. simulans but absent in D. pseudoobscura, and 3 sites with editing specifically present in D. melanogaster. The most parsimonious interpretation is that the brain editome in Drosophila is expanding during evolution (Fig 2). We did not find any convincing case that editing was detected in both D. melanogaster and D. pseudoobscura but was absent in D. simulans, suggesting that the established editing events, at least for the set we examined here, are well maintained by natural selection during evolution. In contrast to previous observations that nonsynonymous editing events were generally non-adaptive in mammals [67,68] and in Drosophila [50], our analysis revealed the nonsynonymous editing events in Drosophila brains were predominantly adaptive. The ratio of nonsynonymous (N) to synonymous (S) editing sites (N/S) in different brain samples of D. melanogaster ranges from 4. 79 with 95% CI (3. 85,6. 20) to 6. 25 (4. 70,8. 67), all of which is significantly higher than the ratio expected under neutrality (3. 80) that was calculated similarly as previously described [67] (Materials and Methods; P < 0. 03 in each library, Fisher’s exact tests; Table 1). In other words, in the brains of D. melanogaster, the rate of nonsynonymous A-to-I editing is significantly higher than the rate of synonymous editing. Given the observed and expected N/S ratios under neutrality (randomness), a conservative estimation is that 20. 7% [with 95% CI (1. 3%, 38. 7%) ] of the N sites in the brains of D. melanogaster might be adaptive (Table 1). Moreover, we obtained higher N/S ratios in each brain library when we increased the cutoff of editing level (S11 Fig). Our analysis is essentially the same as the classical dN/dS (the number of nonsynonymous changes per nonsynonymous site over the number of synonymous changes per synonymous site) test of DNA sequences in molecular evolution [85], and provides compelling evidence that the nonsynonymous editing events in Drosophila brains are overall beneficial and favored by natural selection. We observed significantly negative correlations between the sequencing coverage (C) and editing level (l) in each brain library or when we pooled the library together (P < 10−10 in each case, S12 Fig). These patterns do not necessarily mean that lowly expressed sites have higher editing levels, but rather suggest that the sites of lower sequencing coverage have stronger sampling bias: editing events in such sites are either not detected or detected with over-estimated editing levels. Furthermore, although the S and N editing sites have comparable coverage (Fig 1D), the S sites are generally edited at lower levels compared to the N sites (Fig 1D). Therefore sampling bias would affect S sites more severe than the N sites, which potentially causes over-estimation of the N/S ratios in the above analysis. However, to what extent the N/S ratio is over-estimated due to sampling bias remains unclear. Our “joint probability method” in detecting editing across multiple libraries allows the editing sites with low coverage or with low editing levels in a single library to be efficiently identified with the aid of information from other libraries (such sites are usually filtered out if only based on information of a single library), and the full list of editing sites across libraries enable us to test whether and how our observed N/S ratios are affected by sampling bias. We conducted simulations with two different methods, both of which considered the observed distribution of sequencing coverage and editing levels among sites. In the first method, we focused on all the high-confidence sites present in a library that had sequencing coverage C ≥ Cmin and editing level l ≥ lmin. In each round of simulation, for a site j that had an observed depth Cj, we randomly sampled Cj reads (with replacement) from all the sequenced reads covering that site and calculated the simulated editing level lsj with the obtained reads of the edited allele, we counted the site only if lsj ≥ lmin, and then we pooled all the counted sites together and calculated the N/S ratio for this round of simulation. For each brain library, we tried different combinations of Cmin (ranging from 5 to 20, with Cmin = 20 roughly accounting for 50% of all the editing sites in a library) and lmin (0. 01,0. 02 and 0. 05) values, and performed the simulations for 1,000 replicates. At different lmin cutoffs, both the observed and simulated N/S ratios were generally higher at lower Cmin values; and the median simulated N/S ratio is usually higher than the observed one, but the extent of elevation is very modest, in most cases much smaller than 10% (see Fig 3A and 3B for all the eight brain libraries at lmin of 0. 02, and S13 Fig for other lmin values). Similar results were obtained when we pooled all the libraries together and performed the simulations (S14 Fig). These results suggest that detection bias of editing sites would slightly increase the observed N/S ratio, however, the over-estimation caused by such a bias is very modest compared to the large difference between the observed and the expected N/S ratio under neutral evolution (Fig 3B). Importantly, the observed N/S ratio is significantly higher than the neutral expectation even after we adjusted the bias (S14 Table). In the second method, we pooled all the eight libraries together and randomly sampled (with replacement) a fraction (f) out of the total reads, and after that, we calculated the N/S ratio for the sites that have simulated editing levels lsj ≥ lmin. We tried different combinations of f (from 0. 05 to 1 with a step size of 0. 025) and lmin (0. 01,0. 02 and 0. 05) values, and performed the simulations for 1000 replicates. In agreement with the first method, the simulated N/S ratios were higher at lower depth (smaller f values). For example, when f is set at 0. 05, which is less than half the size of a library we sequenced, the N/S ratio would be elevated by roughly 10% due to sampling bias (Fig 3C). However, the simulated N/S ratios approached to the observed N/S ratio (calculated based on the pooled libraries) rapidly with increasing f (Fig 3C). In summary, the simulations revealed that N/S ratios tend to be overestimated at lower sequencing coverage, however, the degree of over-estimation was usually small given our sequencing depth. Our conclusion that the N/S ratios in the brain libraries are significantly higher than the neutral expectation is not affected by the possible detection bias. It is notable that several X-linked genes such as cac, CG42492 and Sh harbor multiple N editing sites (12,10 and 3 respectively) while very few S editing sites (1,1, and 0 respectively). In fact, the N/S ratio is substantially higher for the X-linked than the autosomal genes in all the eight brain libraries of D. melanogaster, indicating the signal of adaptation is generally stronger for the X-linked genes in Drosophila brains (S15A Fig). This observation is essentially congruent with the fast-X evolution observed for nonsynonymous DNA mutations under positive selection, by which the X-linked advantageous effect is more readily manifested compared to the autosomal counterparts [86,87]. Furthermore, we identified a set of brain-expressed gene candidates that had N editing sites favored by natural selection, which were termed “PSEB” (Positively Selected Editing in Brains) gene set (we required each gene to have a N/S editing ratio > 5; totally 223 genes met this criteria and 49 of them were X-linked; S15 Table). There are 683 high-confidence editing sites in the PSEB genes in brains of D. melanogaster, including 447 N and only 8 S sites, yielding a N/S ratio of 55. 9. 80 (35. 9%) of the PSEB genes are overlapping with the type III genes that have editing events positively selected by Yu et. Al. [50]. The N editing sites in the PSEB genes, which are significantly enriched in chemical and electrical neurotransmission pathways (Fig 4), are very likely the targets of positive selection. Notably, the expression levels of the PSEB genes are significantly higher than the non-PSEB genes (S16 Fig), which suggests that the higher N/S ratios in the PSEB genes are not likely caused by sampling bias of the S sites. Moreover, we ranked all the genes with editing events in each brain library with increasing expression level (RPKM) and equally divided those genes into “Highly” and “Lowly” expressed groups. In both groups, the N/S ratios are significantly higher for the editing sites in the PSEB compared to non-PSEB genes (S16 Table). Furthermore, we conducted random sampling simulations and confirmed that neither the higher N/S ratios in the PSEB genes nor lower N/S ratios in the non-PSEB genes was caused by detection bias (Fig 3D). Notably, the editing levels of the PSEB sites are higher than the remaining sites in all the brain libraries, and four of them are statistically significant (S17 Fig), which further suggests that editing in these sites are favored by natural selection. We observed higher N/S ratios in the sites that have editing events commonly observed across Drosophila species, partially caused by the bias that synonymous adenosine sites are less constrained during evolution. However, after we contrasted the observed N/S ratios to the expected values (4. 18) calculated from the evolutionarily conserved adenosine sites (Materials and Methods), we still detected signals of adaptation in the sites with conserved editing events between D. melanogaster and D. simulans: the N/S ratio ranges from 7. 11 (4. 97,10. 95) to 8. 73 (5. 90,14. 29) across libraries, and all ratios are significantly higher than the neutrally expected ratio (P < 0. 005 in each comparison, Fisher’s exact tests; Table 2). We also observed higher N/S ratios for the conserved editing sites in X chromosome compared to those in the autosomes (S15B Fig). Importantly, the adaptation signals are primarily detected in the PSEB genes, and considerably lower N/S ratios (ranging from 2. 73 to 3. 60) were observed in the non-PSEB genes (S11 Table), suggesting the conserved N sites in the non-PSEB genes are unlikely favored by natural selection. Analogous results were obtained for the conserved editing sites between D. melanogaster and D. pseudoobscura (S12 Table). Altogether, our results suggest that the N editing sites in the PSEB genes are favored and maintained by natural selection while the S editing sites, which are putatively neutral, might degenerate during long-term evolution, which generates even higher N/S ratios in the sites with editing events conserved between species [50,68]. Next we evaluated the effect of local nucleotide contexts on the comparison of observed vs. expected N/S ratios. For each of the 2,114 high-confidence editing sites, we extracted the upstream and downstream 3 nucleotides (we also used other number of nucleotides and obtained similar results), counted the number of nucleotide at each position (S17 Table), and developed a position probability matrix (S18 Table). Consistent with observations in primates that local nucleotide contexts affect the editing efficiencies [60,88], we also found that G immediately upstream a focal editing site was generally not favored, the nucleotide immediately downstream the editing site was slightly biased toward G, and other flanking nucleotides were generally not important, although the overall patterns of preferences were weak in D. melanogaster (Fig 5A; we presented the frequencies of the tri-nucleotides centered with the editing sites and background adenosines in S19 Table). We scanned mRNAs of the edited genes with the position probability matrix and scored each 7-mer sequence that was centered with adenosine (Materials and Methods). We chose the score cutoff that specified the top 90% quantile of the high-confidence editing sites (-0. 622), and found 75. 4% of the background 7-mer sequences were above this score cutoff (Fig 5B). For the background adenosine sites with scores above the cutoff, the expected N/S ratio was 3. 11 (even lower than 3. 80, which was the expected N/S ratio based on all the adenosine sites), significantly lower than the observed N/S ratio (4. 62) for all the high-confidence editing sites that have scores above the threshold (P < 0. 0001, Fisher’s exact test). Therefore, our conclusions based on the comparisons between the observed and expected N/S ratios are not affected after considering effect of the local nucleotides contexts on editing. We calculated the divergence of all the protein-coding genes between D. melanogaster and D. simulans (Materials and Methods), and we found PSEB genes have significantly lower dN and dS values compared to the non-PSEB genes expressed in brains of D. melanogaster (dN is 0. 0107 ± 0. 0012 vs. 0. 0193 ± 0. 0003 for PSEB vs. non-PSEB, and dS is 0. 1025 ± 0. 0025 vs. 0. 1262 ± 0. 0005 for PSEB vs. non-PSEB; P < 10−8 in both comparison, KS test), possibly due to the anti-correlation between expression levels and evolutionary rate [88–91] since the PSEB genes are usually expressed at higher levels (S16 Fig). These observations are consistent with the hypothesis that A-to-I editing increases the mutation sequence space of protein-coding genes. Next we ask whether we can detect higher editing densities in the evolutionarily conserved genes. We ranked all the protein-coding genes expressed in Drosophila brains with increasing dN values between D. melanogaster and D. simulans and grouped the genes into 20 bins. The editing density of the N sites (≥5X coverage) was significantly inversely correlated with the dN value in each bin (rho ranges from -0. 824 to -0. 734, P < 0. 001 in each library, S20 Table; see Fig 6A for 1- to 5-day female and male brains). We observed similar patterns when we ranked all the possible nonsynonymous adenosine sites (≥ 5X coverage was required) with increasing phyloP scores (higher scores mean higher conservation levels) and grouped them into 20 bins: in each library, the density of N sites is significantly positively correlated with the median phyloP score of that bin (rho ranges from 0. 711 to 0. 832; P < 0. 001 in each library, S21 Table; Fig 6B for 1- to 5-day female and male brains). Analogous but weaker correlations were observed for the S sites when we grouped the genes with dN value (rho ranges from -0. 560 to -0. 402 for each individual library, P < 0. 1 in each library, S20 Table) or grouped the synonymous editing sites with phyloP scores (rho ranges from 0. 519 to 0. 738; P < 0. 05 in each library, S21 Table). However, one potential pitfall of our analysis is that the conserved genes (with low dN) or sites (with high phyloP scores) usually have higher expression levels [88–91], which would lead to detection bias as higher editing densities were found in genes (or sites) with higher sequencing coverage (S4A and S4B Fig). To exclude such a possibility, in each library, we ranked all the nonsynonymous adenosine sites with increasing sequencing coverage (≥ 5X was required) and grouped them into 20 bins. Within each bin, we further divided the sites into two equal-sized subgroups based on the phyloP scores. In each library, the editing density in the nonsynonymous adenosine sites is significantly higher in the conserved subgroup compared to the non-conserved subgroup (P < 0. 001 in each comparison, paired t test, Fig 6C and S18 Fig). Analogous patterns were observed for the synonymous editing sites as well (S19 Fig). Therefore, the elevated occurrences of A-to-I editing events in the evolutionarily conserved genes (or sites) in Drosophila are not likely caused by detection bias due to gene expression levels (or bias of sequencing coverage). Altogether, our results further support the adaptation hypothesis of RNA editing [12–14] since the sequence space of the evolutionarily conserved genes is generally inaccessible through DNA mutations, while RNA editing provides an epigenetic approach to expand proteomic diversity temporally and spatially. So far our analysis revealed prevalent beneficial editing sites in brains of Drosophila and the majority of them were enriched in the PSEB genes. Next we asked whether we can observe similar patterns in whole bodies of D. melanogaster adults. We deep sequenced the poly (A) -tailed transcriptomes of female and male adults from five strains of D. melanogaster that were collected from five continents [79] (Materials and Methods). We sequenced 15. 4–28. 6M reads that were mapped on the reference genome in each library (S1 Table), and the median sequencing coverage on an exonic site in a library ranges from 17 to 31 reads in female, and ranges from 11 to 21 in male adults (S20 Fig). We masked all the SNPs in these five and other related strains which were sequenced previously [79], so that we only focused on the DNA sites that were exclusively adenosines across the five strains of D. melanogaster. For each site with A-to-G difference in a strain k, we calculated the probability that this site was edited Pk (E1), and then we calculated the joint probability that this site was edited in at least one strain P (E1). We analyzed female and male adults independently and required each site to have at least 10 raw reads in each library. At FDR of 0. 05, we identified 910 candidate editing sites in female and 1,458 candidates in male adults in exons. We obtained a false positive rate of 4. 17% (26 out of 624) in female and 2. 99% (29 out of 969) in male adults with the w1118 vs. Adar5G1 mutant analysis. After removing the false positive sites in Adar5G1, we obtained 875 exonic editing sites in female (S22 Table) and 1,422 exonic sites in male adults (S23 Table) (719 overlapped between female and male adults). We obtained mixed results in detecting signals of adaptation when we compared the observed N/S ratios to the expected under neutrality (3. 80) in both female and male adults (compared to the neutral expectation, the observed N/S ratio was significantly higher in strain B12 and N10, significantly lower in T07 and ZW155, and no significant difference was observed in I17, Table 4). As observed in brains, the N/S ratio was significantly higher in PSEB genes compared to neutral expectation in all the five strains (P < 0. 001 in each comparison; Table 4). In contrast, in all the five strains the observed N/S ratio was significantly lower in the non-PSEB genes compared to neutral expectation (P < 0. 001 in each comparison; Table 4). The patterns held when we separated the genes with editing events into “Highly” and “Lowly” expressed groups based on their expression levels (S24 Table). We also conducted simulations by randomly sampling the reads covering each site with increasing cutoffs of coverage (Cmin) or editing level (lmin), and our simulation results constantly confirmed these observed patterns (S21 Fig). Our results indicate the N/S ratio is much higher in editing sites in the PSEB genes while strong purifying selection is acting on editing sites in the non-PSEB genes. We uncovered significant but not high correlations in editing levels between strains in female (pairwise Pearson’s r was 0. 718 ±0. 018, S25 Table) and male adults (r was 0. 804 ±0. 013, S26 Table), suggesting that editing levels were variable across strains of D. melanogaster as previously observed [41]. Notably, we found a considerable number of sites at which editing events was readily detected in certain strains while absent in the other strains. We hypothesize that editing in such sites might be polymorphic in the populations, although reliably identifying such sites is challenging. We took a multiple-step procedure to identify the polymorphic editing sites from all the sites we detected (875 in female and 1,422 in male adults). First, we filtered the sites that have editing events detected in D. simulans or D. pseudoobscura. Second, we required a site to have editing reliably detected in at least one strain k with Pk (E1) > 0. 999 and to have no editing detected in at least one strain m. Third, we calculated Pm (D0), the probability that the editing was not detected at depth Cm due to sampling bias or sequencing error in the strain m. Finally we calculated the joint probability P (D0) if no editing was observed in multiple strains at that site (Materials and Methods). For each site, a required parameter in calculating Pm (D0) in the strain m which has no editing detected is the authentic editing level in that strain. We used two levels of stringency to calculate Pm (D0). Level I: at each site, we assumed editing level is identical across all the strains, and the averaged level from the strains with reliably editing signals was used. Level II: we assumed the editing level is 0. 05 in the strain with no editing detected if the averaged level from the detected strains is greater than 0. 05, and the averaged level was used if it was smaller than 0. 05. We conducted the analysis in females and males independently. Finally under Level I we obtained 165 and 179 candidate polymorphic editing sites in females and males respectively (57 sites overlapped). The averaged editing level for a site (based on the strains with reliable editing events) is 0. 125±0. 010 in females and 0. 095±0. 008 in males (we obtained 117 and 125 candidate sites under Level II in females and males respectively, and 44 sites were overlapping; S27 and S28 Table). Interestingly, among these putatively polymorphic editing sites, we did not find signal of adaptation in editing of the PSEB genes (in females the N/S ratio was 4. 67 and 3. 5 under Level I and II respectively; and in males the N/S ratio was 3 and 2. 25 under Level I and II respectively; Table 5). Moreover, we observed even lower N/S ratios in editing of the non-PSEB genes (1. 47 and 1. 48 under Level I and II respectively in females; 1. 96 and 1. 74 under Level I and II respectively in males; Table 5). We observed similar patterns when we individually examined the N and S sites with respect to the number of strains in which such editing events were detected (S29 Table). In contrast, we observed significantly higher N/S ratios in the editing sites of PSEB genes that were fixed in D. melanogaster. To reliably detect the sites with editing events fixed in the populations of D. melanogaster, we first identified the sites at which the probability of editing in each strain Pk (E1) > 0. 95 (k was B12, I17, N10, T07 and ZW155; female and male adults were studied separately). Due to the small number of strains we used, we further sequenced the transcriptomes of female and male adults of D. simulans and required the orthologous sites to be edited in the same gender of D. simulans. We obtained 181 sites (101 in PSEB and 80 in non-PSEB genes, S30 Table) fixed in female and 373 editing sites (194 in PSEB and 179 in non-PSEB genes, S31 Table) fixed in male adults (171 overlapped between female and males). The editing levels are higher in the fixed sites compared to the polymorphic ones: (the averaged editing level per site is 0. 149±0. 014 in female and 0. 197±0. 011 in male adults). In Table 5 we show the fixed editing sites have significantly higher N/S ratios in PSEB genes (44 in females and 34. 2 in males) compared to the neutral backgrounds (4. 18, calculated with the adenosine sites that are in the genes expressed in adults and conserved between D. melanogaster and D. simulans). Furthermore, in the PSEB genes the N/S ratio is significantly higher in the fixed editing sites compared to the sites showing polymorphic editing patterns in females and males after adjusting the background difference (Table 5). Nevertheless, strong signals of purifying selection were observed in the fixed editing sites in the non-PSEB genes (Table 5). We also did not find significant difference in the N/S ratio between the fixed and polymorphic editing sites in the non-PSEB genes in both female and male adults (Table 5). Previous studies have demonstrated that mutations influence mRNA secondary structures and hence the efficiency of RNA editing in natural populations of D. melanogaster [41,92]. Here we ask whether we can find SNPs associated with the variations in editing levels across the five strains of D. melanogaster. To increase the statistical power, we only focused on the 58 editing sites that are polymorphic (under Level I) in both female and males (Materials and Methods). Finally, we found 39 out of the 58 sites are associated with SNPs that are within 10kb flanking the editing sites, and on average each editing site is associated with 22. 4 ± 3. 7 SNPs, and the nearest distance between a SNP and the editing site is 770 ± 140 nts (S32 Table). Meanwhile, we conducted the same analysis on the 171 sites that have editing events fixed in both females and males, however, we did not observe any of those sites have editing levels associated with SNPs under the same criteria. These comparisons suggest cis-regulatory elements affect the levels or status of editing across strains, as previously observed [41,84,92]. Compared to hard-wiring a particular codon change in the genome, A-to-I editing may give the flexibility to quickly respond to environmental stress and adjust the activity of final protein product accordingly [15]. For example, the A-to-I editing level can be regulated by external stimuli [93,94], energy and nutrient [95], and hypoxic conditions [96,97]. Notably, temperature increases would reduce the thermo-stabilities of mRNA secondary structures and down-regulate the expression level of Adar in Drosophila, both of which would reduce the global editing efficiency [16]. The editing levels of Adar and a handful of other genes (totally 54 sites) have been examined under elevated temperatures [16]. Nevertheless, how temperature affects the editing sites at the genome-wide level remains unclear. Herein, we compared the editing levels in the brains of 1- to 5-day-old fly adults constantly raised at 25°C with those in the brains of flies raised at 25°C but treated at 30°C for 48 hours (i. e. , B2 vs. B4 for female, and B6 vs. B8 for male brains of D. melanogaster in Table 1). The changes of editing levels were significantly positively correlated between D. melanogaster and D. simulans (S22 Fig; see Table 2 for information of D. simulans), consistent with previous observation that editing level changes under higher temperature were evolutionarily conserved[16]. The editing levels in general decreased under elevated temperatures for all functional categories of editing events, however, the decrease in editing levels was considerably smaller for the nonsynonymous compared to the silent sites (Fig 7A). Interestingly, when we folded the flanking sequence of each editing event (100 nts at each side), we uncovered the nonsynonymous editing sites were located in more stable secondary structures than all the silent events (P < 10−16, Kolmogorov-Smirnov test, Fig 7B). This comparison suggests that Drosophila has evolved mechanisms to maintain or even enhance the levels of some nonsynonymous editing sites under elevated temperature, which might be important for temperature adaptation since they expanded proteomic diversities. We observed high correlations in editing levels between pairwise brain libraries (Pearson’s r ranges from 0. 849 to 0. 933, S33 Table), and interestingly, the highest correlation coefficients were usually observed between brains of flies that were maintained at the same accommodation conditions but not of the same gender. The pattern was more pronounced when we clustered the samples based on the editing levels of the high-confidence sites (totally 391 sites that had at least 20 raw reads in each brain library): Females and males of 1–14 day old (B1 and B5, see Table 1 for annotations) clustered together; females and males of 1–5 day old that were constantly raised at room temperature (B2 and B6), treated at 30°C for 14 hours (B3 and B7), and treated at 30°C for 48 hours (B4 and B8) always clustered together (Fig 7C; similar patterns were observed when the coverage cutoff was set 15 or 25 raw reads, S23A Fig). Moreover, we did not find any site with editing level significantly different between female and male brains of D. melanogaster under the same accommodation conditions (i. e. , B1 vs. B5, B2 vs. B6, B3 vs. B7 or B4 vs. B8) after multiple testing corrections. Therefore, our results indicate that temperature plays a more important role than gender effect in shaping the global brain editomes. However, we found the species effect is generally stronger than the temperature effect on the editing levels when we clustered the samples of D. melanogaster and D. simulans with 289 sites that have at least 20 raw reads in each library (Fig 7D) or when we clustered the samples of D. melanogaster and D. pseudoobscura with 152 sites that have at least 20 raw reads in each library (S23B Fig): the samples of the same species always clustered first, then the temperature conditions, and the gender effect was still very weak. Gene expression plasticity is a strategy organisms evolved for adapting to new environments [97]. Yet it remains elusive whether (and how) RNA editing and gene expression plasticity coordinately participate in temperature stress responses. We detected hundreds of genes that were significantly differentially expressed in the brains of D. melanogaster that were constantly raised at 25°C compared to those raised at 25°C and treated at 30°C for 48 hours (B2 vs. B4 for female, and B6 vs. B8 for male brains). The down-regulated genes under elevated temperatures were enriched in the “oxidative phosphorylation” pathway in both female and male brains, while the up-regulated genes were enriched in the “ATP binding, ” “translation” and “response to temperature” functional categories in female brains and in the “ATP binding” and “regulation of transcription” pathways in male brains (S34 and S35 Tables), which suggested gene expression plasticity was involved in the temperature stress responses but in a sexually dimorphic manner [98,99]. Interestingly, in both female and male brains of D. melanogaster, the changes in editing levels of the nonsynonymous sites were weakly but significantly positively correlated with the changes in expression levels of the host genes under various cutoffs of expression levels (S24 Fig), suggesting these correlations were not artifacts caused by gene expression cutoffs. For example, a nonsynonymous A-to-I editing site (chrX: 1781840) in Adar mRNA causes a Ser (AGU) to Gly (IGU) change (abbreviated as S>G change), and a previous study [16] and our data both indicated that the editing level of this S>G change (S5A Fig) and expression level (S25 Fig) of Adar mRNA were reduced in both female and male brains of D. melanogaster under elevated temperatures. In contrast, for synonymous sites we did not observe significant correlations between changes in editing levels and gene expression levels (S24 Fig). Altogether, these results suggest that the nonsynonymous editing events might interplay with gene expression changes in temperature adaptation, although detailed mechanisms remain to be further explored.
By extensively characterizing RNA editing sites in brains of D. melanogaster and two related Drosophila species, we identified a considerable number of N sites in Drosophila brains that were adaptive and maintained by natural selection during evolution. Our analysis revealed the N/S ratios in the editomes of Drosophila brains were significantly higher than the neutral expectation. In contrast, we did not observe such a pattern in the editing sites of different developmental stages or whole flies of D. melanogaster that were identified by the modENCODE Project [47] (Fig 8A) or re-analyzed by Ramaswami et al. [52] (S26 Fig). We also obtained mixed results when we compared the overall N/S ratios to the expected ratio under neutral evolution in female and male adults from five strains of D. melanogaster (Table 4). We found the significantly higher N/S ratios in the whole editomes of brains are mainly contributed by the N sites in the PSEB genes, which are favored by natural selection. The N/S ratio for editing sites in the PSEB genes is significantly higher than the neutral expectation in most developmental stages except in early embryos (0–16 hours) or larvae (Fig 8A), while N/S ratio for editing sites in the non-PSEB genes is lower than neutrality in all those samples (Fig 8A). Importantly, in brains, ~60% of the N sites were contributed by PSEB genes, while only ~5% of the S sites were from PSEB genes, which considerably elevated the overall N/S ratio in the brain editomes (Fig 8A). In contrast, in the late embryo, pupae and adults, although the N/S ratios for the PSEB sites were significantly higher than neutral expectation, less than 40% of the N sites were contributed by PSEB genes in each stage/tissue (Fig 8B), and hence, the signatures of adaptation in PSEB sites are masked by the non-PSEB sites when pooling all the editing sites together (Fig 8A). These patterns were constantly observed when we independently considered the “Highly” and “Lowly” expressed genes that have editing events in these samples (S27 Fig). We also observed a trend that the editing levels in the PSEB sites were increased during Drosophila development (S28 Fig). Furthermore, we observed a significant positive correlation between the expression level of Adar and the number of editing sites during Drosophila development in the modENCODE Project [47] (Spearman’s rho = 0. 413, P = 0. 023; S29A Fig). And the expression level of Adar was also positively correlated with the cumulative editing levels (i. e. , summing up the editing levels of all the sites, as following Ref. [63]) in the modENCODE samples (rho = 0. 487, P = 0. 006; S29B Fig). Importantly, the expression level of Adar is higher in brains compared to all the modENCODE developmental stages (Fig 8C), and accordingly, we observed the highest number of editing sites in brains compared to other developmental stages (S30 Fig). Notably, in early embryos (0–16 hours) and larvae, either Adar was lowly expressed, or only a few PSEB genes were expressed in these stages (Fig 8C), which putatively explains why we did not detect significantly higher N/S ratio for editing sites in PSEB genes in these stages (Fig 8A). Taken together, our results demonstrate that the expression level of Adar, together with the expression profiles of the PSEB genes that have editing sites favored by natural selection, are important in shaping the overall N/S ratios in the global editomes at different developmental stages (or tissues) of D. melanogaster. We were able to predict secondary structures for ~74% of the exonic editing sites, including ~22% as long-range pseudoknots, which were often very challenging to find. Our results help understanding the molecular basis by which the editomes are regulated and maintained. First, our analysis reveals a considerable number of editing sites could be clustered due to promiscuous editing of multiple adenosines by ADAR simultaneously. Second, we found a significantly higher proportion of the N sites were located in stable hairpin structures of mRNAs than the silent sites (55. 0% vs. 41. 7%, P < 0. 001, Fisher’s exact test), and similarly, the flanking sequences of the N sites have significantly lower MFE compared to those of the silent sites (Fig 7B). Since ADAR recognize double-stranded RNAs to exert editing, these findings provide the structural basis for the observed higher editing levels and excessive occurrences of the N sites. The density of editing events was significantly higher in the evolutionarily conserved genes than in the non-conserved ones (Fig 6C, S18 Fig), which supports the hypothesis that editing increases the mutation sequence space. However, another competing explanation for this observation is that the secondary structures (stable hairpins or pseudo-knots) of mRNAs, which are ADAR substrates for RNA editing, constrain the evolutionary rates of these genes. Indeed, when we separately calculated the evolutionary rates between D. melanogaster and D. simulans for the CDS regions that formed secondary structures and harbored editing events in Drosophila brains (termed “structured”) and the remaining CDS regions after masking the secondary structures (“structure masked”), we uncovered both dN and dS were significantly lower in the “structured” than “non-structured” CDS regions (P < 10−10 for both dN and dS comparisons, S31 Fig; totally 172 genes were included in the comparisons). Nevertheless, the length of the “structured” CDS regions in general only account for 4. 36 ± 0. 24% of the full lengths of CDS regions, and we still observed significantly lower dN and dS values in the “structure masked” CDS regions compared to the genes without editing events (“Unedited”) in Drosophila brains (P < 10−10 for both dN and dS comparisons, S31 Fig). Hence, although the editing-associated secondary structures considerably constrain mRNA sequences, they generally have negligible impact on the evolutionary rates of the total CDS regions. In summary, our results support the hypothesis that A-to-I editing increases the proteomic diversity for the genes that are highly conserved due to functional constraints. It is worth noting that compared to humans [67,68] and macaques [59,60], we observed strong signals of adaptation in the editing sites in Drosophila, especially in the brain editomes. Furthermore, the A-to-I editing events we identified in Drosophila were significantly enriched in evolutionarily conserved genes while the editing sites in human coding regions showed an opposite pattern [67]. The different observations between Drosophila and primates might be shaped by the difference in the underlying molecular mechanisms and selective forces. First, there are two catalytically active ADAR enzymes (ADAR1 and ADAR2) in primates and both enzymes are expressed in many tissues, which potentially cause promiscuous editing events that might be neutral or deleterious [59,60]. Second, the targets of editing are mainly repetitive non-coding sequences in primates [13]. Although the coding editing events conserved between human and mouse are adaptive [68], the majority of the editing events in coding regions might be solely by-products of the over-activity of ADARs and hence selected against [67]. In contrast, there is only one Adar locus in Drosophila [37,38], which is predominately expressed in the nervous system and preferentially edits pre-mRNAs of neural genes [39]. Third, the effective population size is much larger for D. melanogaster than primates, which makes natural selection more efficient in the former than in the latter [100]. Therefore, the adaptive editing events, once originated, will be more effectively spread and fixed in Drosophila than in primates. On the other hand, the detrimental effects incurred by RNA editing, will be more efficiently selected against in D. melanogaster than in primates. Besides providing proteomic diversity, our results also suggest that mRNA editing interplayed with gene expression plasticity to fine-tune gene expression activity under temperature stress responses, which supports previous hypothesis that RNA editing might be a driving force for environmental adaptation [15,16]. Interestingly, we also found several editing events in D. melanogaster that compensated for the G-to-A DNA mutation in the D. melanogaster lineage after splitting with its sibling species (S32 Fig), which suggests that RNA editing events are advantageous because they reverse the deleterious effects caused by G-to-A DNA mutations as previously proposed [67,101]. Taken together, our evolutionary analyses, combined with our functional genomic studies, shed new light on the molecular mechanisms and functional consequences of RNA editing in Drosophila.
Flies were grown in 12 hour light: 12 hour dark cycles at 25°C. The ISO-1 strain of D. melanogaster, the sim4 strain of D. simulans, and one lab strain of D. pseudoobscura were gifts from Dr. Andrew G. Clark’s lab at Cornell University. In the temperature stress experiments, 1- to 5-day-old flies were transferred from 25°C incubators to 30°C incubators and treated for 14 hours or 48 hours, and the humidity and light conditions were maintained at the same levels. The 1- to 5-day-old and 1- to 14-day-old flies were separately sexed, and the brains were dissected in RNAlater solution (Ambion). We also separated the heads and bodies of the female or male adults of the sim4 strain of D. simulans with a fine sieve and extracted the total RNAs from heads and bodies of each gender. Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer' s protocol. Poly (A) + mRNA was isolated from 15 μg total RNA with oligodT25 DynaBeads (Thermo Fisher). Next, the mRNA was fragmented and size selected from 40 nts to 80 nts by 15% TBU gels. Following 3' dephosphorylation, 3' ligation with a 3' adaptor, 5' phosphorylation and 5' ligation with a 5' adaptor, and size-selected mRNA fragments were reverse transcribed with SuperScript III (Invitrogen). The sequence of the 5' adaptor was 5' GUUCAGAGUUCUACAGUCCGACGAUC3' and the 3' adaptor was 5' TGGAATTCTCGGGTGCCAAGG3' . All cDNA was amplified by 14 PCR cycles with Phusion High-Fidelity DNA polymerase (NEB) with the TruSeq index adapters, and the products within the correct size ranges were collected from 20% TBE gels for the quality tests (Fragment Analyzer, Agilent Technologies) and sequencing (Platform: Illumina HiSeq 2500; read length: 50 bp, single-end). The 3' adaptor sequences were clipped by the Cutadapt program [102]. The remaining reads were aligned to the reference genome of D. melanogaster (r6. 04), D. simulans (r1. 4) or D. pseudoobscura (r3. 2) using STAR [103] (mapping statistics were summarized in S1 Table). The genome sequences and annotations of the three Drosophila species were downloaded from FlyBase (www. flybase. org). First, for each brain library, we employed the GATK RNA-Seq variant calling pipeline [78] to detect the A-to-I editing events, where the transcriptomic base is Adenosine and the sequencing read variant base is Guanine (Inosine). We were able to identify 1,531 editing sites at this preliminary stage. To retrieve the editing sites that were potentially excluded by the pretty strict GATK pipeline, we pooled all the editing sites in four previous studies [47–50] and the GATK candidates altogether, constructing a list of 5,925 candidate sites. Second, for each candidate site in a brain library k, we discarded the reads with mapping quality lower than 10 and the reads with mismatches other than A-to-G, and extracted the sequencing coverage (Ck) and the number of edited allele that shows A-to-G difference (Lk) with SAMtools [104] and calculate the probability that a site is edited in a library k given the observed sequencing data Pk (E1) = 1 − Pk (E0), where Pk (E0) is the probability the A-to-G difference is solely caused by sequencing error with a rate of ε. We define Pk (E0) =∑i=LkCk (Cki) ∙εi∙ (1−ε) Ck−i. The Illumina HiSeq platform generally has an error rate ε0 ranging from 0. 2% to 0. 6% [105–107] and we used ε0 = 0. 5% in the analysis. As we focus on the sites with solely A-to-G but not A-to-C or A-to-T mismatches, the error rate of A-to-G would be scaled to ε=ε0/31− (2/3) ε0≈0. 00167. The probability that a site is edited in none of the n libraries is thus P (E0) =∏k=1nPk (E0), where n is the number of brain libraries, and accordingly, the probability that this site is edited in at least one of these libraries is P (E1) = 1 − P (E0). For each of the 5,261 sites that have sequencing coverage in our brain libraries of D. melanogaster, we calculated P (E1) and P (E0), and corrected for multiple testing with Benjamini & Hochberg method [108]. The procedures are summarized in Fig 1A. The functional annotations of the editing sites were conducted with the software SnpEff [109]. The transcriptome data for heads of the Adar5G1 mutant and paired wild type w1118 of D. melanogaster were downloaded from NCBI SRA under accession numbers SRR629970 and SRR629969 [51]. We mapped the transcriptome data to the reference genome of D. melanogaster using STAR [103] and extracted the A and G alleles with SAMtools [104]. To estimate the false positive rates of the editing sites we detected, we first examined the number of the editing sites we identified in D. melanogaster that also had editing events detected in the heads of w1118 strain (N1), and then we counted the number of these shared sites that also showed A-to-G difference in the Adar5G1 mutant (N2). The false positive rate is estimated by N2/N1 as conducted previously [50,51]. We folded all the transcripts in genes expressed in brains of D. melanogaster with RNALfold [110] and identified the exonic editing sites that were located in the stable local hairpin structures (z score < -1. 5, ΔG < -15 kcal/mol, and the stem length > 50 nts). To obtain the expected numbers of editing sites in the hairpin structures by randomness, we randomly sampled the same number of exonic editing sites under study with replacement for 1000 replicates and counted the median and 2. 5% and 97. 5% quantiles in the simulations. For each introinc editing site, we folded the flanking sequences for each site (100 nts at each side) with RNALfold and examined whether this intronic site was located in the stable hairpin structures with the same criteria. For each of the editing events detected in D. melanogaster, the flanking sequences (100 nts at each side) were folded In the long-range pseudoknot searches, the full-length pre-mRNA or flanking sequences (2,000 nts at each flank, totally 4,000 nts for long pre-mRNAs) were folded with RNAfold [110] and putative long-range pseudoknots were parsed if the pairing region in each stem of the pseudoknot was > 40 nts and the distance between the two pairing stems was > 100 nts. We also fold the flanking sequences for each editing site (100 nts at each side) and calculated the MFE (kcal/mol) with RNAfold [110]. The protein and CDS sequences of D. melanogaster (r6. 04), D. simulans (r1. 4) and D. pseudoobscura (r3. 2) were downloaded from FlyBase. The reciprocal best orthologous genes were obtained based on pairwise BLASTP [111] between D. melanogaster and D. simulans (r1. 4), and between D. melanogaster and D. pseudoobscura. The protein sequences of the orthologous genes were aligned with the clustalw [112] program, and the CDS alignments were produced with the tranalign [113] program based on the corresponding protein alignments. The yn00 program in the PAML [114] package was employed to calculate the dN and dS values for each gene between D. melanogaster and D. simulans. The phyloP score for each site of D. melanogaster were downloaded from UCSC Genome Browser (genome. ucsc. edu). For each editing site in D. melanogaster, we employed two complementary approaches to search for the orthologous sites in D. simulans and D. pseudoobscura. First, we used liftOver [83] to convert the genomic coordinates of the orthologous sites in coding and non-coding regions between D. melanogaster and D. simulans, or between D. melanogaster and D. pseudoobscura as previously conducted [51] (termed “g_align” approach). The pairwise genome alignments between D. melanogaster and D. simulans, and between D. melanogaster and D. pseudoobscura were downloaded from UCSC Genome Browser (genome. ucsc. edu) and used to identify the evolutionarily conserved adenosine sites. Second, we parsed the genomic co-ordinates with the pairwise CDS alignments between D. melanogaster and D. simulans, and between D. melanogaster and D. pseudoobscura (termed “c_align” approach). The g_align approach efficiently identified both the coding and non-coding orthologous sites, and the c_align approach is powerful in identifying orthologous sites in coding regions between distantly-related species. Among the 2,114 high-confidence editing sites in brains of D. melanogaster, we identified 1,499 orthologous sites in D. simulans that were also adenosines and had sequencing coverage in at least one brain library of D. simulans (1,443 by g_align, 577 by c_align, and 521 by both), and 892 sites in D. pseudoobscura that were also adenosines and had sequencing coverage in at least one brain library of D. pseudoobscura (707 by g_align, 527 by c_align, and 342 by both). To exclude SNPs in the RNA editing characterization, we also deep sequenced the genomic DNA of the same strains of D. simulans (the median coverage per site is 46) and D. pseudoobscura (the median coverage per site is 47) that were used for RNA-editing detection. We mapped the DNA reads on the reference genomes with BWA[115], excluded reads with mapping quality lower than 10, and called the SNPs with SAMtools (313,133 SNPs in D. simulans and 489,828 SNPs in D. pseudoobscura). After masking the SNPs, for each site in each species, we calculated P (E1), the joint probability that this site is edited in brains. At FDR of 0. 05, we identified 996 sites edited in D. simulans (947 by g_align, 495 by c_align and 446 by both), and 451 sites edited in D. pseudoobscura (340 by g_align, 326 by c_align, and 215 by both). Our experimental designs (brain samples of the same gender and the same age under the same accommodation conditions for different species) and the combinations of g_align and c_align approaches enabled us to identify more evolutionarily conserved editing events compared to Yu et al. [50] which mainly focused on the conserved editing events in the coding regions. For an editing site with coverage Cm and editing level lm in a library m, no edited (G) allele would be detected if the edited RNA molecules were not sampled or all the edited signals were abolished by sequencing errors. Therefore, we estimated the probability of observing zero edited reads at this site in a sample m as Pm (D0) =∑i=0Cm (Cmi) ∙lmi∙ (1−lm) Cm−i∙εi∙ (1−ε) Cm−i, where ε is the scaled sequencing error rate (0. 00167). Then the joint probability that this site is edited but not detected in n libraries would be P (D0) =∏m=1nPm (D0). The expected N/S ratio for the editing events under neutral evolution was calculated with a similar procedure described previously [67]. Briefly, for each adenosine site in the CDS regions of the genes that have at least one editing event in brains of D. melanogaster, we tested whether it cause an amino acid change (nonsynonymous, N) or not (synonymous, S) when edited. To obtain the expected N/S ratios for the editing events that are evolutionarily conserved, we conducted similar analysis only on the conserved adenosine sites between D. melanogaster and D. simulans, or between D. melanogaster and D. pseudoobscura, as previously conducted between human and mouse [68]. The N/S ratio expected under neutral evolution is 3. 80 for the editing sites in D. melanogaster, 4. 18 for the editing sites with events conserved between D. melanogaster and D. simulans, and 7. 11 for the editing sites with events conserved between D. melanogaster and D. pseudoobscura. We obtained the 95% confidence intervals by random sampling the same number of observed editing sites with replacement and calculated N/S ratios from the simulated data for 1000 replicates. The detailed procedures of the two methods in evaluating the effect of detection bias on N/S ratio are fully described in TEXT. The processes of these simulations (and other simulations and statistical tests in this study) were performed using R (www. r-project. org). For each of the 2,114 high-confidence editing sites, we extracted the upstream and downstream 3 nucleotides flanking this site, counted the number of nucleotide at each position of the 7-mer, and developed a position probability matrix. Sequence logo for this motif was generated with WebLogo (http: //weblogo. berkeley. edu). In order to get a well-controlled set of genomic background sites, we scanned the mRNA regions for other adenosine sites in genes where the 2,114 editing events were detected. Then we scored both the 2,114 high-confidence editing sites and the background adenosine sites with the same position probability matrix. For a given sequence N-3N-2N-1AN1N2N3 (Ni ∈ (A, T, C, G) ), the score for this sequence was calculated as ∑ilog2Pi (Ni) 0. 25, where Pi (Ni) is the probability of observing base Ni at position i (-3, -2, -1,1, 2,3) based on the position probability matrix. We chose the score cutoff that specified the bottom 10% quantile of the high-confidence editing sites (-0. 622). With that score cutoff, we chose the sites with scores above this cutoff (-0. 622) in the high-confidence editing sites and background adenosine sites, and calculated the N/S ratio in each dataset. The 1–14 day old female and male adults from five strains (B12, I17, N10, T07 and ZW155) of D. melanogaster were sexed, the poly-A tailed mRNAs were isolated from females and males independently with the procedures described above. The libraries were prepared with NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina and the sequencing was carried out at Illumina HiSeq 2500 (read length: 100 bp, paired-end). These strains, kindly provided by Dr. Andrew G. Clark at Cornell University, were originally collected from five continents [79]: Beijing, China (abbreviated B); Ithaca, NY USA (abbreviated I); Netherlands, Europe (abbreviated N); Tasmania, Australia (abbreviated T); and Zimbabwe, Africa (abbreviated Z). These flies were grown in 12 hour light: 12 hour dark cycles at 25°C. The A-to-G editing sites in the exonic regions were called with the “joint probability” method as above described, and the SNPs in these five strains and other 79 related strains based on the whole-genome re-sequencing information [79] were masked in the down-stream analysis. Finally we identified 875 editing sites in female and 1,422 exonic sites in male adults. To identify the editing sites with putative polymorphic events across these five strains, we first filtered the sites that have editing events detected in brains of D. simulans or D. pseudoobscura, or whole files of D. simulans. And then we required a site to have editing reliably detected in at least one strain k with Pk (E1) > 0. 999 and to have no editing detected in at least one strain. Next we calculated Pm (D0), the probability that the editing was not detected at depth Cm due to sampling bias or sequencing error in a strain m as above described. Finally, we calculated the joint probability P (D0) if no editing was observed in multiple strains at that site. Under Level I in calculating Pm (D0), we assumed the strain m that has no editing event detected have the same editing level at that site as in the other strains which had reliably editing event detected (the mean value was used). Under Level II, we used 0. 05 as the editing level in calculating Pm (D0) if the mean editing level from the other strain is >0. 05. To detect the editing events that were fixed across these five strains, we employed two complementary approaches: First we identified the sites at which the probability of editing in each strain Pk (E1) > 0. 95 (k was B12, I17, N10, T07 and ZW155; we studied female and male adults separately); and then we sequenced the transcriptomes of female and male adults of D. simulans (the sim4 strain) and required the orthologous sites to be edited in the same gender of D. simulans. To identify the SNPs that are associated with the variation of the editing levels in female or male adults across these strains, we only focused on the editing sites that are polymorphic in both female and males (totally 58 sites) or fixed in both females and males (171 sites). For each site in females, we retrieved the SNPs within 10kb regions that flanked the focal editing sites in the five strains and conducted the association test between the editing level in each strain and the genotypes of a SNP (we assumed the reference allele as 0 and the alternative allele as 1). We also performed the same analysis in male adults. We required the reference and alternative allele of a SNP to be associated with editing levels of the same site in the same direction in females and males and P < 0. 05 in both tests. Total RNA from the female or male brains was prepared independently using Sanger sequencing procedures. The total RNA was treated with RNase-free DNase I (Invitrogen) to remove genomic DNA. Reverse transcription was performed using random primers, and cDNA was amplified using target-specific primers (the primers sequences are presented in S36 Table). The final PCR products were sequenced with the Sanger method at the Ruibiotech Sequencing Company. qRT-PCR was performed with SYBR Green Master Mix (Thermo Fisher) in a 20 μL reaction volume and monitored on a StepOnePlus Real-Time PCR System (Thermo Fisher). rp49 was used as the internal control. The primers for the real-time PCR assay are listed in S37 Table. The raw NGS reads for each gene were counted with the htseq-count program [116], gene expression levels were normalized, and differentially expressed genes were detected with the edgeR package [117]. The RPKM for each gene or half-gene was calculated with CuffLinks [118]. All of the gene ontology analyses were performed by DAVID [119] and all the brain-expressed genes were used as background list. The sequence data in this study have been submitted to the NCBI Sequence Read Archive (SRA, http: //www. ncbi. nlm. nih. gov/sra) under accession number SRP074828 and SRP068882. All other relevant data are within the paper and SI files.
Adenosine-to-inosine (A-to-I) RNA editing is an evolutionarily conserved mechanism that alters RNA sequences at the co-transcriptional or post-transcriptional level. RNA editing is hypothesized to facilitate adaptation in that it expands the transcriptomic and proteomic diversity. However, evidence for adaptation of RNA editing at the whole editome level is still lacking. In this study we systematically identified A-to-I RNA editing sites in female and male brains of three Drosophila species at different temperatures. With evolutionary analysis from different perspectives, we provide lines of evidence to demonstrate that the nonsynonymous editing sites in Drosophila brains are generally adaptive. The signals of adaptation for the editing sites are significantly enriched in genes related to chemical and electrical neurotransmission. We show that the RNA editing events might interplay with gene expression plasticity in temperature stress responses. Furthermore, we demonstrated that the expression level of Adar, together with the expression profiles of a set of genes that have editing sites favored by natural selection, were important in shaping the overall selective patterns of the global editomes at different developmental stages (or tissues) of D. melanogaster. Altogether our results support the hypothesis that A-to-I editing provides a driving force for adaptive evolution in Drosophila from different aspects.
invertebrates
glycosylamines
computational biology
animals
invertebrate genomics
animal models
drosophila melanogaster
model organisms
experimental organism systems
genome analysis
drosophila
adenosine
research and analysis methods
genomic libraries
gene expression
evolutionary genetics
insects
animal genomics
arthropoda
biochemistry
rna
rna editing
nucleic acids
nucleosides
genetics
biology and life sciences
genomics
evolutionary biology
glycobiology
organisms
Adaptation of A-to-I RNA editing in Drosophila