Introduction

RNAs undergo many post-transcriptional changes (McCown et al., 2020), with RNA editing being one of the more widely observed mRNA modifications. RNA editing alters single nucleotides in precursor mRNA and is most often catalyzed by the ADAR (adenosine deaminase acting on RNA) family of enzymes that bind double stranded mRNA and irreversibly deaminate specific adenosines into inosines (A-to-I editing) (Keegan et al., 2001; Nishikura, 2010; Rieder and Reenan, 2012; Shevchenko and Morris, 2018). Inosines are recognized as guanines by the translation machinery, and edits within mRNA coding (CDS) exons can alter the protein’s amino acid sequence. In addition, editing of 5’ or 3’ untranslated (UTR) exons can regulate mRNA splicing, translation, stability or microRNA binding (Berkovits and Mayr, 2015; Borchert et al., 2009; Lev-Maor et al., 2007; Lopdell et al., 2019; Mazloomian and Meyer, 2015; Nakano et al., 2017; Rueter et al., 1999; Shevchenko and Morris, 2018; Stellos et al., 2016; Yuan et al., 2023). The levels of RNA editing (the fraction of total mRNA in an individual cell with a specific edit) depend on editing site, developmental stage and cell type (Birk et al., 2023; Buchumenski et al., 2017; Graveley et al., 2011; Kliuchnikova et al., 2020; Lundin et al., 2020; Ma et al., 2002; Picardi et al., 2017; Sapiro et al., 2019; Wahlstedt et al., 2009). Although editing occurs across the animal kingdom, the biological significance of most editing events is unclear. Mammalian genomes encode two ADAR proteins (ADAR1 and ADAR2), with ADAR1 implicated in innate immunity (Keegan et al., 2023). The majority (>95%) of editing found in humans is mediated by ADAR1 and occurs within copies of the primate-specific non-coding Alu mobile retroelements that comprise ∼11% of the human genome (Athanasiadis et al., 2004; Bazak et al., 2014; Eisenberg and Levanon, 2018; Levanon et al., 2004). ADAR2 has broader functional impacts through its ability to recode amino acids within a small number of proteins (Keegan et al., 2023; Pinto et al., 2014). RNA editing is particularly robust in neurons and several targets have been shown to regulate membrane excitability and synaptic function (Alon et al., 2015; Brija et al., 2023; Brusa et al., 1995; Hoopengardner et al., 2003; Jepson and Reenan, 2009; Palladino et al., 2000a; Yablonovitch et al., 2017).

In contrast to mammals, the Drosophila genome encodes a single ADAR that is functionally more similar to mammalian ADAR2 (Palladino et al., 2000b). Thousands of RNA editing sites have been identified in Drosophila using pooled RNA samples from whole animal or brain homogenates, enriched GFP-labeled cell populations obtained by flow cytometry, or through computational predications of intron-exon base pairing (Graveley et al., 2011; Hoopengardner et al., 2003; Keegan et al., 2023; Kliuchnikova et al., 2020; Mazloomian and Meyer, 2015; Ramaswami et al., 2013; Rodriguez et al., 2012; Sapiro et al., 2019; Stapleton et al., 2002; St Laurent et al., 2013). Mutations in the sole Adar gene in Drosophila cause behavioral dysfunction, temperature-dependent seizures and neurodegeneration in surviving adults (Jepson and Reenan, 2009; Keegan et al., 2011; Khan et al., 2020; Li et al., 2014; Palladino et al., 2000a; Robinson et al., 2016). Although ADAR is important in the Drosophila nervous system, the critical editing targets and functional impact of most edits remain unclear.

To examine the landscape of RNA editing more systematically, we determined the editing profile for Drosophila larval motoneurons (MNs) that are commonly used to study neuronal development and glutamatergic synapse biology (Harris and Littleton, 2015). Given the robust imaging and electrophysiology toolkits available to examine synaptic function at larval neuromuscular junctions (NMJs), combined with GAL4 drivers for specific MN subtypes (Aponte-Santiago et al., 2020; Pérez-Moreno and O’Kane, 2019), these neurons provide an excellent model to define the functional impacts of RNA editing on neuronal communication. Two well-characterized glutamatergic MN subtypes innervate most larval bodywall muscles. Tonic-like Ib neurons target single muscles and form weaker presynaptic active zones (AZs) that display synaptic facilitation, while phasic-like Is neurons synapse onto multiple muscles and form stronger AZs that undergo synaptic depression (Genç and Davis, 2019; He et al., 2023; Jetti et al., 2023; Kurdyak et al., 1994; Lnenicka and Keshishian, 2000; Medeiros et al., 2023; Newman et al., 2017; Sauvola et al., 2021). These MN subtypes also have diverse biophysical and morphological features (Aponte-Santiago and Littleton, 2020), providing an opportunity to identify the major targets for RNA editing and any cell-type specific edits. We recently determined the transcriptomes for larval Ib and Is MN subtypes by performing single cell Patch-seq RNA profiling from hundreds of individual Ib and Is neurons (Jetti et al., 2023). This allowed identification of differentially expressed genes (DEGs) that contribute to the unique properties of these neurons. In the current study, we extend this analysis to define the RNA editing landscape and single cell editing rules for these two neuronal populations. These data identify the major targets for RNA editing in the two cell types and reveal a highly stochastic RNA editing process that occurs across individual MNs from both populations.

Results

Single motoneuron RNA editing suggests functional diversification of select number of target neuronal proteins

Canonical RNA editing occurs when ADAR deaminates adenosines, converting them into inosines that are interpreted as guanosine by the translation machinery (Figure 1A). To determine specific sites that are edited in Drosophila 3^rd instar larval motoneuron (MN) subpopulations, we analyzed single cell RNA sequencing data obtained by Patch-seq profiling from 105 individual MN1-Ib (hereafter referred to as Ib) and 101 individual MNISN-Is (hereafter referred to as Is) neurons labeled with GFP using GAL4 drivers specific to each cell type (Jetti et al., 2023). Candidate RNA editing sites were obtained using GATK (McKenna et al., 2010; Poplin et al., 2017) and SAMtools (Danecek et al., 2021) to compare genomic DNA sequenced from parental strains to their single cell RNA transcriptomes to identify base pair mismatches that did not represent genomic SNPs (Figure 1B, C). To focus on high-confidence RNA editing events and further exclude potential sequencing errors that could appear as very low abundance edits, DNA-RNA mismatches were filtered and excluded if the edited site was observed at low levels (<10% editing) or present in only a few neurons (<10 Ib or Is cells). After filtering, 316 high confidence canonical A-to-I editing sites were identified across both MN transcriptomes (Figure 1D, Supplemental Tables 1, 2). A number of C-to-T and G-to-A transversions were also observed (Figure 1D, Supplemental Tables 3, 4), suggesting noncanonical editing occurs in MNs as seen in other Drosophila cell types (Sapiro et al., 2019; Stapleton et al., 2006). Although noncanonical editing has been described (Bhakta and Tsukahara, 2022; Gott and Emeson, 2000; McCown et al., 2020; Pecori et al., 2022), the biological significance and molecular pathways mediating these edits are less clear. We largely focused on canonical A-to-I editing given it is likely to have a greater impact on MN function, in contrast to noncanonical editing that appears less impactful based on the target sites identified (see analysis of noncanonical edits below).

Figure 1.

The 316 canonical edits were distributed across 210 genes and included previously unannotated sites based on Flybase (Supplemental Table 2). 55% (175/316) of canonical edits occurred in mRNA coding sequences (CDS), while 36% (115/316) were in 3’ UTR and 8% (26/316) in 5’UTR (Figure 1E), consistent with prior editing distributions in Drosophila (St Laurent et al., 2013). The editing rate within each region of the mRNA was normally distributed (p < 0.05, Kolmogorov-Smirnov test). Mean editing rate in the CDS (0.71) was significantly higher (ANOVA F statistic = 20.171, degrees of freedom (df) = 2, Tukey p < 7.8×10^-10) than 3’UTR (0.59) and was also higher than 5’UTR (0.64) editing (Figure 1F). For CDS edits, 58% (102/175) resulted in a missense amino acid substitution. Given edits that result in significant amino acid changes are more likely to alter protein function, we analyzed these edits in more detail. Missense changes were classified as significant if they altered the residue classification between (but not within) positive (R, K) or negative (D, E), hydrophobic (A, V, I, L, M, F, Y, W), polar (N, Q), potential phosphorylation sites (S, T), or structural resides (C, G, P, H). Based on these criteria, 59% (60/102) of missense edits resulted in significant alterations to the affected amino acid. The most common missense change was L>M, while the most common disruptive alteration was change of an existing residue into a glycine (S>G, E>G, R>G, and D>G, Figure 1G). No canonical edit was identified that induced a premature stop codon. However, an edit in Cnep1r2, which encodes a serine/threonine phosphatase involved in nuclear pore formation and lipid metabolism (Jacquemyn et al., 2021), changed a stop codon into a tryptophan residue. The unedited transcript encodes a protein that is conserved with its mammalian homologs through the C-terminus, while the edited transcript generates a modified protein with 14 additional C-terminal amino acids.

Editing rules from single neuron analysis

Given the RNA editing data derive from single neurons with potentially distinct RNA editing patterns, we were interested in defining general rules for single cell editing compared to those observed using bulk RNAseq approaches. Early studies of editing rates for mammalian GluR2 AMPA receptors observed either a lack of editing or complete editing of GluR2 mRNA depending on developmental stage, suggesting an all-or-none editing model. However, recent analyses indicate editing rates for most sites is more variable (Graveley et al., 2011; Rodriguez et al., 2012; Sapiro et al., 2019; St Laurent et al., 2013). To examine the editing profile across every editing site in individual Ib and Is neurons, the average fraction of edited reads in single cells was determined (Figure 1F). Editing rates spanned from the 10% threshold to 100% across both neuronal populations. The average editing rate per site was 66% and was not significantly different between Ib and Is neurons (Figure 1F, p = 0.94, Student’s t-test, t statistic = 0.8, df=630). Consistent with prior studies in adult Drosophila CNS (Hoopengardner et al., 2003; Palladino et al., 2000a; Sapiro et al., 2019), mRNAs encoding regulators of neuronal excitability and function were over-represented as editing targets in larval MN transcriptomes (Figure 2A). The largest number of edits (15 sites) within a single gene occurred in the mRNA encoding Complexin (Cpx), a key regulator of synaptic vesicle (SV) fusion (Buhl et al., 2013; Cho et al., 2015, 2014; Huntwork and Littleton, 2007; Jorquera et al., 2012). Other highly edited mRNAs included those encoding voltage-gated ion channels (Cac, Para, Shab, Sh, Slo), ligand-gated ion channels (nAChRalpha5, nAChRalpha6, nAChRalpha7, nAChRbeta1, Rdl) and regulators of synaptic transmission (EndoA, Syx1A, Rim, RBP, Vap33, Lap).

Figure 2.

RNA editing events that induce functional changes to a transcript would be predicted to be expressed in most cells. Although no edit site was detected in 100% of sampled neurons, often due to low read coverage at specific edit sites in any given cell, some edits were found in as many as 87% (175/206) of neurons. On average, an edited site was detected in ∼19% (∼40/206) of neurons (Figure 2B). Edit sites observed in the largest fraction of sampled neurons are shown in Figure 2C. Edits in the 5’ and 3’UTR were more common in these cases, with edits in ATPsynCF6 and Syx1a found in >80% of individual neurons. There was also a tendency for edit events to occur in proximity to other edits. Out of 316 sites, 53% had a nearby edit within 100 nucleotides and 26% had a nearby edit with 10 nucleotides (Figure 2D), suggesting favorable ADAR binding to some dsRNA regions increases editing rates at nearby adenosines (Buhl et al., 2013). Indeed, 12 edits with an average editing rate of 45% occurred in a stretch of 118 nucleotides in the 3’UTR of Cpx (3R:4301140-4301258). In summary, these data indicate RNA editing in larval MNs acts on a small number of target mRNAs and ranges from low rates of editing to sites where the mRNA is fully edited.

Identification of editing sites that alter conserved amino acids

Amino acid changes that occur in conserved regions of a protein are more likely to have significant impact on function. We filtered significant missense changes caused by editing and focused on those occurring in conserved residues or protein domains based on homology across other invertebrate or vertebrate species. 43 editing-induced significant missense changes were identified in conserved areas across several protein classes, including those regulating membrane excitability, synaptic transmission, and cell signaling (Table 1). 49% (21/43) of the edited residues altered amino acids that were conserved in vertebrates and invertebrates, while the remaining 51% (22/43) changed residues conserved across other insect homologs. Notably, many of the edits occurred within annotated protein domains. Given these edits alter conserved amino acids and are likely to have a larger impact on larval MN function, they were examined in more detail.

Table 1.

Multiple subunits of the nicotinic acetylcholine receptor (nAChR) family, which function as ligand-gated cation channels mediating fast excitatory neurotransmission, were found to undergo editing. One edit in nAChRalpha5 converts a threonine to alanine (T791A) in the 4^th transmembrane domain and is edited at very high levels (Ib: 99%, Is: 96%, Figure 3A), suggesting the edited version represents the majority of nAChRalpha5 protein in MNs. Edits in nAChRalpha6 (Figure 3B) and nAChRbeta1 (Figure 3C) altered conserved residues in the extracellular ligand-binding domains of the proteins. Editing also occurred in several proteins regulating neuronal excitability, including E-to-G and Y-to-C edits in the ATPα sodium pump (Figure 3D), a highly edited S-to-G site (Ib: 86%, Is: 91%) in the 9^th transmembrane segment of the voltage-gated calcium channel Cacophony (Figure 3E), and multiple edits in the ion transport domain of other voltage-gated ion channels, including Shab, Sh, Slo and Para (Figure 3F).

Figure 3.

Beyond neuronal excitability, mRNAs encoding proteins regulating the SV cycle were also editing targets (Table 1). Three edit sites within the C-terminal domain of the Cpx7A splice isoform were identified and previously described (Brija et al., 2023). Editing of adjacent adenosines encoding the N130 amino acid produce N130S, N130G or N130D variants that alter Cpx’s ability to regulate SV fusion (Brija et al., 2023). Editing was also observed in a key regulator of endocytosis, Endophilin A (EndoA), resulting in a change from a positively charged K to a negatively charged E residue at amino acid 140 within the F-BAR lipid-binding domain. Given the K residue is conserved in all Endophilin homologs, one would predict this charge substitution alters EndoA function in endocytosis within MNs. Editing of another SV endocytosis regulator, the Clathrin AP180 adaptor homolog Lap, occurred at high rates (Ib: 88%, Is: 95%) and resulted in a T-to-A change in a conserved residue. Edits were also identified in the AZ proteins Rim and RBP that altered C-terminal residues conserved in other insect homologs. Similarly, an edit in the VhaAC45RP subunit of the SV proton pump resulted in neutralization of a negatively charged E residue conserved in other insect homologs.

Within the group of edits predicted to significantly impact protein function, 67% (29/43) were previously annotated on FlyBase as editing sites, while 33% (14/43) were not previously described. Newly identified sites could be enriched targets for 3^rd instar larval MNs or missed in prior work. Among the new edits was a missense D-to-G change in Bsk (Figure 3G), a conserved serine/threonine protein kinase central to the stress-induced c-Jun N-terminal kinase (JNK) pathway. The edited D residue is universally conserved in JNK kinase homologs and resides within the core of the kinase domain. Editing from D to G occurs at moderate levels (Ib: 34%, Is: 30%) and is predicted to be disruptive to the protein’s function. Another novel editing site was identified in the spliceosomal protein Beag (Figure 3H). Beag forms a complex with Smu1 and additional alternative splicing proteins and is required for larval synapse development and function (Beck et al., 2012). The edited serine residue is conserved across insect species and is edited at high rates within MNs (Ib: 72%, Is: 57%). The S-to-G substitution occurs in the RED-like protein C-terminal domain, a region involved in nuclear localization and binding to other nuclear-localized proteins (Ulrich et al., 2016). A third newly identified edit site results in an E-to-G substitution in the highly conserved endoribonuclease RNase Z (Figure 3I) involved in tRNA processing (Saoura et al., 2017). The E744 residue is universally conserved in invertebrate and vertebrate homologs and localizes downstream of its metallo-beta-lactamase C-terminal domain. These data indicate RNA editing alters highly conserved residues across several proteins, including those involved in membrane excitability and synaptic transmission that are likely to impact the development or function of larval MNs.

Characterization of highly edited transcripts and editing stochasticity

We next examined sites with a very high editing rate that also induced missense amino acid substitutions to identify edits that may cause larger functional impacts on MNs similar to mammalian AMPA receptor editing. 27 sites in 18 genes were found to be edited at >90% levels in Ib or Is neurons that underwent editing at the location (Supplemental Table 5). 85% (23/27) of these sites were previously annotated on FlyBase, consistent with these targets being more commonly identified given their high editing rate. Highly edited sites were overrepresented in genes regulating membrane excitability and synaptic function. 10 sites in 9 genes displayed 100% editing in at least one of the two neuronal populations, including: (1) quiver, a trafficking regulator of K⁺ and nACH receptors; (2, 3) nAChRalpha5 and nAChRalpha6, subunits of the nicotinic acetylcholine receptor; (4) Prosap, the Drosophila SHANK synaptic scaffolding homolog; (5) Rdl, a GABA receptor; (6) Nckx30C, a Na⁺/K⁺-dependent Ca²⁺ exchanger; (7) Shab, a delayed rectifier K⁺ channel; (8) Ythdc1, a regulator of alternative splicing; and (9) RNaseZ, a pre-tRNA endoribonuclease. Although these sites were highly edited, most sites showed more variable patterns of editing, with some neurons editing the site at high levels while others had lower rates or no editing. Adar mRNA levels were not significantly different between Ib and Is neurons (Figure 4A, p = 0.55, adjusted p-value from RNA sequencing Wald statistic (Jetti et al., 2023)). Similar to prior studies (Graveley et al., 2011; Lundin et al., 2020; Picardi et al., 2017; Wahlstedt et al., 2009), Adar mRNA expression level at the resolution of single cells did not correlate with overall editing rates (black line in Figure 4B). Indeed, neurons with low or high Adar transcript levels were able to edit specific sites up to 100%. However, we noticed that edit sites in abundantly expressed mRNAs (expression level in the top 30% as noted by red dots in Figure 4B) generally had lower editing rates. Given that 4000-fold differences in gene expression can occur across transcripts within larval MNs (Jetti et al., 2023), we hypothesized that abundant transcripts might result in lower editing if ADAR activity becomes rate-limiting. Indeed, mRNA expression level was inversely correlated to the fraction of reads edited at individual sites for both Ib and Is subtypes (Ib: Pearson’s r = −0.59, r² = 0.35; Is: r = −0.58, r² = 0.34, Figure 4C). We conclude that editing levels are not directly limited by ADAR expression but rather by target gene expression, suggesting ADAR can become rate-limiting when high levels of an editable pre-mRNA transcript are present.

Figure 4.

Although Adar mRNA levels did not correlate with editing rate, it is possible the expression of other RNA binding proteins or editing factors impact the rate of editing. If so, individual neurons that express higher levels of these factors would be expected to display higher overall editing rates. To examine this possibility, editing rate for every target site was plotted for each individual Ib and Is MN (Figure 4D). Neurons from both populations displayed a range of average editing rates across all target sites that varied from ∼30% to >80%, indicating the presence of “low” and “high” editors (Figure 4D). Although it is unclear what drives these differences in editing rate for neurons of the same cell type, future analysis of DEGs between these groups might identify candidate proteins that impact editing across an otherwise homogenous neuronal population. Another possibility for differential editing is that cells with the highest editing rate edit fewer transcripts overall but at a higher rate per transcript. Consistent with this model, the mean fraction of reads edited versus the number of edits across individual neurons displayed a negative correlation (Ib: Pearson’s r = −0.33, r² = 0.11; Is: r = −0.50, r² = 0.25, Figure 4E). We hypothesized that an edit found in most cells would potentially have a high editing rate. However, a negative correlation was observed, suggesting the more cells an edit was found in, the lower the mean editing rate (Figure 4F). Beyond variability in overall editing rate, some individual edit sites displayed a highly stochastic cell-to-cell editing rate whereby some cells expressed a mixture of edited and unedited mRNA, whereas others had all or none editing for that same target site. Several representative regions showing stochastic editing rates are shown for three such edits in Cpx (Figure 4G). In summary, editing rates for most target sites are variable at the level of single neurons and individual edited adenosine residues, with only a small number of sites undergoing complete editing.

Comparison of RNA editing between Ib and Is motoneuron populations

Given structural and functional differences between tonic Ib and phasic Is MNs, we hypothesized unique RNA editing targets may help establish their distinct features. However, the overall landscape of RNA editing across Ib and Is neurons was similar at the level of edit site location (5’UTR, CDS, 3’UTR) and editing rate (Figure 1F). In addition, the number of edit sites and fraction of reads edited per cell showed similar profiles (Figure 5A). A comparison of editing rate at individual sites indicated most were edited to the same level between Ib and Is neurons (Pearson’s r = 0.73, r² = 0.53, Figure 5B). The fraction of Ib and Is MNs that expressed individual edits was also correlated (Pearson’s r = 0.83, r² = 0.69), although several edits were only detected in one of the two cell types (red and blue dots in Figure 5C). Though editing rules were generally similar between the two populations, several differences emerged. The lower tail of the distribution for the fraction of edited reads displayed a sharper drop across Ib cells (Figure 4D), indicating a subset of these neurons have distinctly lower editing rates than observed in Is neurons. Similarly, the upper tail of the distribution revealed a subset of Is neurons with higher editing rates than observed within the Ib population (Figure 4D).

Figure 5.

To identify sites with significantly different editing rates between subtypes, differentially expressed mRNAs were removed from the analysis. Following DEG filtering, 26 RNA edit sites across 22 genes were identified that displayed statistically significant editing differences (Figure 5D, E, Supplemental Table 6). Although 3’UTR edits comprised only 36% of all editing sites (Figure 1B), 46% (12/26) of differentially edited sites between Ib and Is neurons were in 3’UTRs. Given 3’UTR editing can modify microRNA binding sites or protein translation rates, it is interesting to consider if post-transcriptional regulation through RNA editing in this region differentially impacts protein production. Of particular interest are two 3’UTR sites in complexin (3R:4301150 and 3R:4301158) that are edited at greater rates in Is neurons, which are known to express lower Cpx protein levels at synapses (Newman et al., 2022). We conclude that differential editing between Ib and Is neurons that alters amino acid identity is rare and unlikely to play a major role in the morphological and functional diversity between these neuronal populations. Differential editing within the 3’UTR of several genes could impact expression of their encoded proteins, but neuronal subtype diversity of Ib and Is cells appears to be driven mostly by their unique DEGs rather than differential RNA editing.

Larval muscles have fewer editing targets and reduced editing rate per site compared to motoneurons

To compare the editing landscape of MNs with a non-neuronal cell type, Patch-seq data generated from larval abdominal muscles 1 and 4 (Jetti et al., 2023), which are postsynaptic targets for Ib and Is MNs (Figure 1C), was analyzed using the same pipeline as described above. Drosophila bodywall muscles are generated from fusion of individual myocytes, forming polynuclear cells containing ∼8-20 nuclei that share a common cytoplasm (Windner et al., 2019). As such, Patch-seq from larval muscles is likely to sample a larger population of editing possibilities compared to single nuclei MNs. Despite being polynucleated, only 37 canonical editing sites were detected in at least 10 muscle samples that had >10% editing rate (Supplemental Tables 1, 2, and 7, Figure 6A). Seven of these sites were uniquely edited in muscles and the remaining 30 were also present in Ib or Is MNs (Figure 6A). Muscle cells displayed similar or lower editing rates than Ib or Is neurons at shared edit sites (Figure 6B, Supplemental Table 2). Muscle-specific editing sites (Supplemental Table 7) were enriched for highly expressed muscle genes, including Act5C (Actin homolog), wupA (Troponin 1 homolog), TpnC47D (Troponin C homolog) and mthl3 (Secretin-like GPCR). Like MNs, editing rate in muscles was variable for specific target sites, ranging from 18% at the lower end to 75% edited (Supplemental Table 7). Compared to an average editing rate per site of 66% for MN targets, the average editing rate for muscle-only targets was 34%. These data indicate larval MNs have more editing sites and higher editing rates than their postsynaptic muscle targets, consistent with RNA editing having a greater impact on neuronal function.

Figure 6.

Noncanonical RNA editing is observed in larval motoneurons but results in fewer significant amino acid changes

Although canonical A-to-I ADAR-dependent deamination is well-known, other forms of mRNA editing have been described and are emerging in importance due to their application in targeted genome editing (Christofi and Zaravinos, 2019; Xu et al., 2022; Zhang et al., 2024). By comparing MN transcriptomes to parental genome sequences, we identified noncanonical ribonucleoside substitutions (non-A-to-I) across several mRNAs (Figure 1D, Supplemental Tables 3, 4). A similar percent of edits occurred in the CDS for canonical (55%) and noncanonical (48%) editing (Figure 7A, B, Figure 1E, F). However, the fraction of silent changes was much higher for noncanonical than canonical edits (70% versus 42%, Figure 7C). Although fewer of these edits resulted in an amino acid change, they generated a distinct pattern of substitutions compared to ADAR-dependent editing (Figure 7 – figure supplement 1A). Unlike canonical editing, noncanonical edit sites were not enriched for neuron-specific genes (Supplemental Table 4, Figure 7 – figure supplement 1B). However, the overall fraction of cells displaying noncanonical RNA edits (Figure 7 – figure supplement 1C) and the likelihood of edits to cluster near each other (Figure 7 – figure supplement 1D) were similar to ADAR-dependent editing. In addition, there was evidence for a rate-limiting noncanonical editing machinery for abundantly edited transcripts (Figure 7D) and a reduced number of overall edits for cells with a high editing rate (Figure 7E). Noncanonical editing also exhibited variability in rate per site across individual cells for both Ib and Is subtypes, with some neurons displaying higher or lower editing rates (Figure 7F). Overall, Ib and Is neurons showed similar numbers of edited sites (Figure 7 – figure supplement 1E), editing rate per individual site (Figure 7 – figure supplement 1F), and fraction of cells containing a specific noncanonical edit (Figure 7 – figure supplement 1G).

Figure 7.

Only a few noncanonical edits were identified that altered conserved residues in key neuronal proteins (Supplemental Table 4), suggesting a limited impact on MN function. One edit of interest was identified in the activity-regulated Arc1 mRNA that encodes a retroviral-like Gag protein involved in synaptic plasticity (Ashley et al., 2018; Hantak et al., 2021). This C-to-T edit resulted in a P124L amino acid change, with ∼50% of Arc1 mRNA edited at this site for both Ib and Is neurons. Arc1 assembles into higher order viral-like capsids (Erlendsson et al., 2020) and P124 resides in the hinge region between the central capsomere structure (amino acids 42-120) and the C terminal end (amino acids 125-203) that protrudes from the core capsomere. P124 is highly conserved in other insect Arc homologs, suggesting the edit could alter capsid formation or structure. We also identified this Arc1 C-to-T edit in publicly available RNAseq datasets from adult Drosophila antennal olfactory neurons (Harrell et al., 2021), with an edit frequency ranging from 20-80% (Figure 7 – figure supplement 2A-C). Although Arc1 is an activity-regulated gene, an increase in neuronal activity driven by overexpression of the inward rectifier Kir2.1 channel did not impact its editing (Figure 7 – figure supplement 2C, p = 0.3792, Student’s t-test, t statistic = 0.8, df = 22). The Arc1 C-to-T edit was also identified in another publicly available RNAseq dataset from adult Drosophila brains (Bukhari et al., 2024). While the edit was observed at ∼50% in control brains, it was absent in a matched Alzheimer’s Disease model expressing pathogenic Tau proteins (Figure 7 – figure supplement 2D), suggesting pathogenic Tau expression alters Arc1 editing. Given Arc1 mRNA levels are upregulated following Tau expression (Bukhari et al., 2024), the higher abundance of the mRNA may become rate-limiting for editing and partially account for the reduction. We also identified a noncanonical G-to-A edit that resulted in a G2E missense mutation in Frequenin 1 (Frq1), a multi-EF hand containing Ca²⁺ binding protein that regulates synaptic transmission (Dason et al., 2009). This glycine residue is conserved in all Frq homologs and is the site of N-myristoylation, a co-translational lipid modification that tethers Frq proteins to the plasma membrane (Farazi et al., 2001). Some Frq family members cycle on and off membranes when this domain is sequestered in the absence of Ca²⁺ (Burgoyne et al., 2004). Although the G2E edit is found at lower levels (20% in Ib MNs and 18% in Is MNs), it could create a pool of soluble Frq1 that alters downstream signaling. Given the overall lack of neuronally-enriched target sites and a larger propensity for generating silent coding changes, we conclude that the widespread effects of canonical A-to-I editing within Drosophila larval MNs likely play a more significant role in generating functional protein diversity than noncanonical editing.

Comparing Patch-seq target sites and editing rate to other larval and adult neuron RNAseq datasets

We next identified RNA editing targets in several other publicly available RNAseq datasets that used pooled larval or adult neurons to compare to those found using our Patch-seq approach. Raw RNAseq datasets were analyzed from three sources: (1) 3^rd instar larval MNs that were sorted by FACS (fluorescence-activated cell sorting) using the type I MN-specific OK6-GAL4 driver (Cypranowska et al., 2024), hereafter referred to as “Larva”; (2) nuclei from adult head neurons FACS sorted using the pan-neuronal elav-GAL4 driver (Sapiro et al., 2019), hereafter referred to as “Adult1”; and (3) total nuclei isolated from adult heads comparing nascent versus polyA+ RNA (Rodriguez et al., 2012); hereafter referred to as “Adult2”. Additional SNP filtering was done by comparing Patch-seq data with genomic sequence from control yw adults from the Adult2 study (Rodriguez et al., 2012). Sites displaying any reads with the edited base were eliminated from further analysis. After this additional SNP filtering, 127 total canonical CDS, 111 canonical UTR and 403 noncanonical CDS sites from the Patch-seq data were examined further (Supplemental Tables 8 and 9).

Similar to Patch-seq editing, we required that >10% of the mRNA was edited at an individual site to be considered a positive hit across the three additional datasets. Among 127 Patch-seq canonical CDS edits, 67 were previously annotated in Flybase (Figure 8A). As expected, all of the previously annotated sites were identified in at least one other RNAseq dataset, with most being detected in all 3 groups (Figure 8B). Among the 60 unannotated sites found using Patch-seq, 58.3% were validated in at least 1 other dataset (Figure 8C, D). The canonical R-to-G edit in the Rim AZ protein was robustly detected in all three groups and represents an example of an unannotated edit supported in all RNAseq datasets examined. From the three datasets, the largest proportion of unannotated edits were in the Larva RNAseq data, suggesting they may be MN specific or developmentally regulated (Figure 8D). These edits include a D-to-G edit in CG42613, a CUB-domain containing membrane protein homologous to the Tolloid-like protein 1 (TLL1) protease, a D-to-G edit in Ankyrin 2 (Ank2), and an E-to-G edit in the Hormone receptor-like 78 (Hr78) protein. Other missense-causing edits in this group included an I-to-V edit in BetaTub97EF and a M-to-V edit in the Shab K⁺ channel. Unannotated edits in the Larva dataset that were present generally had higher editing rates (>50%) in the Patch-seq data, suggesting some unannotated canonical CDS edits observed only in with Patch-seq might be specific to the MN1-Ib and MN-Is neuron subtypes and masked in pooled samples. Several canonical CDS edits from the Patch-seq data were not identified in the other datasets, including the S-to-G edit in Beag, the D-to-G edit in Bsk, and the E-to-G edit in RNaseZ. These edits could be neuron subtype-specific or represent false positives. Further studies will be required to confirm the validity of this smaller subset.

Figure 8.

To determine if highly edited sites identified via Patch-seq displayed similarly high rates in other datasets, the editing frequencies from pooled 1b and 1s MNs from this study were compared to compiled frequencies from the other three RNAseq groups. For previously annotated sites, a positive correlation in editing frequency was observed across all three datasets (Larva vs. Patch-seq annotated: Pearson’s r = 0.78, r² = 0.61; Adult1 vs. Patch-seq: Pearson’s r = 0.70, r² = 0.49; Adult2 vs. Patch-seq: Pearson’s r = 0.67, r² = 0.45, Figure 8E-G, blue dots). For example, the edit in nAChRalpha5 resulting in a T791A change was highly edited across all datasets (>90%), while the ATPα sodium pump edit converting E-to-G had lower editing levels (5-20%), consistent with Patch-seq. The correlation plots also identified several sites whose editing frequencies increased in adult brain samples. The R to G edit in nAchRbeta1 and the N-to-S edit in Para displayed higher editing levels in the Adult1 (nAChR: 63%, para: 80%) and Adult2 (nAChR: 66%, para: 82%) datasets relative to Patchseq (nAChR: 36%, para: 68%) and Larva (nAchR: 16%, para: 48%) data (Figure 8E-G). A weaker positive correlation was observed for unannotated sites (Larva vs. Patch-seq unannotated: Pearson’s r = 0.40, r² = 0.16; Adult1 vs. Patch-seq: Pearson’s r = 0.47, r² = 0.22; Adult2 vs. Patch-seq: Pearson’s r = 0.27, r² = 0.07, Figure 8E-G, orange dots). This observation suggests some unannotated sites are likely Ib and Is MN-specific and the editing rate becomes diluted in pooled RNAseq datasets, while others may represent false positives. The Patch-seq editing frequencies for unannotated sites (0.74) and annotated sites (0.71) confirmed in at least one dataset were significantly higher compared to unannotated edits identified only with Patch-seq (0.55, ANOVA F statistic = 5.732, df = 2, Tukey p = 0.0042), indicating the higher the editing rate found in Patch-seq, the greater the likelihood the edit was detected in other published studies (Figure 8H).

Among the 111 Patch-seq canonical UTR edits, 16 were previously annotated in Flybase and 95 were unannotated (Figure 8 – figure supplement 1A). As with CDS sites, 100% of the previously annotated sites were validated, with the majority being detected in all 3 groups (Figure 8 – figure supplement 1B). Among the 95 unannotated sites, 63.2% were validated in at least 1 dataset (Figure 8 – figure supplement 1C). Of these confirmed unannotated edits, 14.7% were detected only in the Larva dataset (Figure 8 – figure supplement 1D). Three Syx 5’ UTR edits belonged to this category, indicating these represent larval and/or MN-specific edits. 11 of the 12 unannotated Cpx 3’ UTR edits were also validated across the other datasets. In terms of UTR editing frequency, there was a positive correlation between Patch-seq and other RNAseq data for previously annotated sites (Larva vs. Patch-seq annotated sites: Pearson’s r = 0.90, r² = 0.81; Adult1 vs. Patch-seq annotated sites: Pearson’s r = 0.90, r² = 0.81; Adult2 vs. Patch-seq annotated sites: Pearson’s r = 0.84, r² = 0.71, Figure 8 – figure supplement 1E-G, blue dots), indicating the relative editing levels at these sites are robust and static across development and neuronal cell types. For example, the previously annotated Slo 3’ UTR 3R:24707423 site was highly edited (80-96%) across all 3 RNAseq groups and correlated with Patch-seq levels (94%). Similarly, the previously annotated EndoA 3’ UTR 3R:18910153 site was edited at moderate levels across the three groups (44-57%) and the Patch-seq data (66%). Similar to CDS edits, editing levels for unannotated UTR sites showed weaker correlation (Larva vs. Patch-seq unannotated sites: Pearson’s r = 0.22, r² = 0.047; Adult1 vs. Patch-seq unannotated sites: Pearson’s r = 0.15, r² = 0.02; Adult2 vs. Patch-seq unannotated sites: Pearson’s r = 0.09, r² = 0.008, Figure 8 – figure supplement 1E-G, orange dots). From this group, several sites displayed consistent editing frequencies, including an edit in the Sh 3’ UTR (30-50% across the three groups and 67% in Patch-seq) and an edit in the nAchRbeta1 3’ UTR (30-70% across the three groups and 50% in Patch-seq). Unlike CDS sites, the Patch-seq editing frequencies for unannotated UTR sites (0.57) and annotated sites (0.51) confirmed by at least one dataset were similar to sites only identified in Patch-seq (0.57, ANOVA F statistic = 0.4950, df = 2, Tukey p = 0.6110, Figure 8 – figure supplement 1H), indicating UTR edits are generally not unique to Ib and Is MN cell types.

All 403 noncanonical edits identified in Patch-seq were previously unannotated, as most editing studies have excluded non-A-to-I editing from analysis. Among the noncanonical edits, 45% (182/403) were confirmed in at least one of the RNAseq datasets (Figure 8 – figure supplement 2A). The largest proportion of confirmed edits (16.4%) were verified in the Adult2 group, with 13.7% of confirmed edits verified only in the Larva dataset (Figure 8 – figure supplement 2B). Noncanonical edits that appeared only in the Patch-seq and Larva datasets included an E-to-Q edit in the RNA binding protein Staufen, a G-to-E edit in the mitochondrial ribosomal protein mRpS5, and an A-to-S edit in the ribosomal protein RpS15. Noncanonical edits in this category that are predicted to be less disruptive included A-to-V and E-to-D edits in the ribosomal protein RpS15Ab and a M-to-L edit in Syx1A. Editing frequencies for noncanonical CDS events between Patch-seq and other RNAseq datasets were poorly correlated (Larva vs. Patch-seq: Pearson’s r = 0.25, r² = 0.065; Adult1 vs. Patch-seq: Pearson’s r = 0.21, r² = 0.045; Adult2 vs. Patch-seq: Pearson’s r = 0.11, r² = 0.012, Figure 8 – figure supplement 2C-E). Similar to canonical CDS edits, noncanonical edits that were more highly edited (>50%) in Patch-seq data were more likely to be detected in other datasets. Indeed, Patch-seq editing frequencies of noncanonical sites confirmed in other datasets were significantly higher (0.68) than sites identified only in Patch-seq (0.51), indicating the higher the editing in larval MNs, the more likely the edit was detected in other datasets (Figure 8 – figure supplement 2F, p < 0.0001, Student’s t-test, t statistic = 7.304, df=368). Several noncanonical editing events also clustered at 50% or 100% editing in the Larva and Adult1 datasets, raising the possibility they correspond to previously unknown SNPs from the parental genotypes used in these studies. Indeed, this clustering was not observed in the Adult2 dataset where parental genomic DNA was sequenced and all SNPs removed. We conclude that although noncanonical editing events are found in all three datasets, the lower correlation with the Patch-seq data suggests some of these sites may represent unknown SNPs or occur at higher frequency in larval MNs.

To examine if editing identified by Patch-seq was disrupted in Adar mutants, editing frequency was compared between adult head samples from control and Adar null mutant datasets (Rodriguez et al., 2012). The vast majority of canonical CDS and UTR editing events identified, including previously unannotated sites, were abolished in Adar mutants (Figure 8 – figure supplement 3A, B, Supplemental Table 8), indicating these are ADAR-dependent. Noncanonical CDS sites were largely unaffected by loss of ADAR as expected (Figure 8 – figure supplement 3C, Supplemental Table 9). To examine if edits identified by Patch-seq were likely to occur co-transcriptionally, editing frequencies between nascent RNA (isolated from adult head nuclei) and polyA+ RNA (isolated from adult heads) were compared using datasets from Rodriguez et al. (2012). Editing levels for canonical CDS and UTR edits were well correlated between nascent and polyA+ mRNA (Canonical CDS: Pearson’s r = 0.90, r² = 0.80; Canonical UTR: Pearson’s r = 0.88, r² = 0.78, Figure 8 – figure supplement 3D-E), suggesting most edits occur co-transcriptionally, consistent with previous findings (Rodriguez et al., 2012). In contrast, noncanonical editing sites were poorly correlated between nascent and polyA+ mRNA (Pearson’s r = 0.34, r² = 0.12, Figure 8 – figure supplement 3F). A subset of these sites had higher editing in polyA+ mRNA (i.e. R-to-H in CG42518, 48% edited in polyA+ versus 0% in nascent RNA), suggesting some noncanonical editing may preferentially occur in processed mRNA following nuclear export. Indeed, proteins implicated in C-to-U and G-to-A editing localize to both the nucleus and cytoplasm in mammalian cells (Anant et al., 2001; Ikeda et al., 2011).

Finally, we examined if any RNA editing events identified by Patch-seq were activity-regulated by taking advantage of Larva RNAseq data obtained following reductions in neurotransmitter release using knockdown of the AZ proteins RBP and Unc13 (Cypranowska et al., 2024). Disruptions to neurotransmitter release under these conditions results in a homeostatic response in which MNs increase their excitability to compensate. Given many editing targets in MNs regulate neuronal excitability or synaptic transmission, changes in editing at these sites might alter MN excitability. Editing frequencies were compared for canonical CDS, canonical UTR and noncanonical CDS sites in larval MN RNAseq datasets from controls or following RBP or Unc13 RNAi (Figure 8 – figure supplement 3G-L). Overall, editing frequency was largely unchanged following presynaptic disruptions to neurotransmitter release (RBP vs. control canonical CDS: Pearson’s r = 0.94, r² = 0.88; RBP vs. control canonical UTR: Pearson’s r = 0.96, r² = 0.91; RBP vs. control noncanonical CDS: Pearson’s r = 0.96, r² = 0.93; Unc13 vs. control canonical CDS: Pearson’s r = 0.94, r² = 0.87; Unc13 vs. control canonical UTR: Pearson’s r = 0.97, r² = 0.94; Unc13 vs. control noncanonical CDS: Pearson’s r = 0.96, r² = 0.93, Figure 8 – figure supplement 3G-L). However, a few sites displayed evidence for a change in editing rate. Among canonical CDS sites, the Y-to-C edit in the ATPα sodium pump was markedly reduced for both RNAi conditions (RBP: 6%, Unc13: 9%) compared to controls (44%). Editing rates for the Shab S-to-G and Cac S-to-G sites were also reduced following Unc13 knockdown. A modest increase in the Rim R-to-G edit was observed in both RBP and Unc13 RNAi knockdown datasets. For UTR editing, a slight increase in the nAchRbeta1 3’ UTR site and a modest decrease in the Slo 3’ UTR editing frequencies were observed in both RNAi datasets. A 3’ UTR edit in Sh was slightly increased following Rbp RNAi. Several noncanonical edits were also affected by presynaptic weakening (Figure 8– figure supplement 3I, L ), including a decrease in an M-to-L edit in AnxB10 in both knockdown conditions. Given Rim, Cac, Shab, Sh, and Slo mRNAs were identified as downregulated DEGs under these hypoactivity conditions (Cypranowska et al., 2024), it will be informative to examine if these changes in target editing frequency contribute to their reduced mRNA expression.

Discussion

Multiple types of ribonucleoside modifications occur across mRNAs, tRNAs, rRNAs and ncRNAs (McCown et al., 2020). RNA editing modifies single nucleotides in mRNA transcripts in a cell-type-specific manner and is one of the more widely observed RNA modifications across metazoans. In this study, we characterized the RNA editing landscape of individual Drosophila 3^rd instar larval MN populations using Patch-seq RNA extraction to define rules for single cell editing and identify the most relevant editing targets. Tonic-like Ib and phasic-like Is MN form well-studied glutamatergic synapses with cell-type-specific electrophysiological and morphological features, providing a model system to investigate how RNA editing contributes to their morphology and function, as well as any role in cell-type-specific diversity. The majority of A-to-I edits identified in our study occurred within mRNA coding regions, with 60 edits predicted to cause significant alterations to protein function. As observed in prior studies (Hoopengardner et al., 2003; Jepson and Reenan, 2009; Palladino et al., 2000a; Rueter et al., 1999; Sapiro et al., 2019), mRNAs encoding ion channels and synaptic proteins were among the most highly edited targets, suggesting editing of these gene classes are likely to fine tune membrane excitability and synaptic output in larval MNs. The overall RNA editing landscape was similar across Ib and Is populations, indicating differences in editing targets or editing rate is unlikely to be a major driver of cell-type diversity for these neuronal subtypes.

Using single cell sequencing from stereotyped and identified single Ib and Is MNs, this analysis revealed a wide variation in average RNA editing levels for specific sites and across individual cells. Indeed, many sites underwent highly stochastic editing when compared across the neuronal population. For example, some neurons displayed ∼100% editing at certain sites, while others displayed no editing for the same target. Such dramatic differences in editing rate at specific target sites is likely to contribute to the heterogeneous features observed within the same neuronal population. At the single neuron level, each cell would be predicted to express a population of both edited and unedited proteins. What drives the variability in RNA editing across individual MNs is unclear. Given ADAR is an enzyme, its affinity for specific dsRNA structures plays a role in establishing editing levels for any specific pre-mRNA target site.

Although Adar mRNA level in individual neurons was not correlated with editing rate, the expression level of the target mRNA undergoing editing was a major factor. Indeed, mRNA expression level was anti-correlated with editing percent, suggesting Drosophila ADAR activity becomes rate-limiting for abundant pre-mRNAs as proposed in mammals (Lopdell et al., 2019). Beyond mRNA expression level, variable expression of other trans-acting factors that regulate ADAR function or pre-mRNA accessibility could modulate editing rate or target site selection across individual neurons as previously hypothesized (Ingleby et al., 2009). Given single cell DNA-RNA mismatches were excluded if fewer than 10 total RNAseq reads covered the site, editing calls for poorly expressed mRNAs within these neuronal populations are not included in this study. We also identified edits that were only present in larval MNs (this study and Cypranowska et al., 2024) and not observed in adult RNAseq datasets, consistent with cell type or developmental stage-specific editing. In terms of activity-regulation, a few targets showed changed in editing frequency when presynaptic neurotransmitter release was reduced.

Although this study defines the RNA editing landscape for Drosophila larval MNs, future experiments will be required to elucidate how specific RNA editing events regulate target gene expression or protein function. Of particular interest will be edits occurring at high levels, such as the 3’UTR edit in ATPsynCF6 and 5’UTR edits in Syx1A, as these may represent edits required for normal transcript localization, expression, or function. For editing events that lead to amino acid change, these substitutions could in principle improve, disrupt, or modify protein function. We failed to observe cases of re-coding an amino acid from a non-conserved residue into a residue conserved in orthologous proteins, suggesting most coding edits with a functional effect are likely to diversify or disrupt protein activity as previously suggested based on comparative genomics (Duan et al., 2023; Liu et al., 2024). The evolutionarily advantage, if any, that RNA editing has in modifying protein function is still unclear. One hypothesis is that RNA editing plays similar roles to alternative splicing. For example, changes in the ratio of Unc13A and Unc13B splice isoforms during Drosophila AZ development regulate the probability of SV release by altering the position of docked SVs near Ca²⁺ channels (Pooryasin et al., 2021; Reddy-Alla et al., 2017). In a similar manner, ratiometric pools of edited versus unedited proteins could provide a mechanism to fine-tune specific neuronal features, particularly for proteins whose abundance is rate-limiting or that act within multimeric complexes. Such a mechanism has been observed for distinct Cpx edit variants (Brija et al., 2023), where individual Cpx proteins are required to clamp multiple assembling SNAREpins during SV priming. In addition, functional properties of the Shaker voltage-gated K⁺ channel are shaped by the combination of individual edits across 4 distinct sites (Ingleby et al., 2009).

Canonical and noncanonical (non-A-to-I) editing was observed in larval MNs, though ADAR-dependent editing resulted in far more edits to neuronal-specific genes predicted to impact protein function. The majority of noncanonical CDS edits caused silent changes, suggesting evolutionary selection against this mechanism as a pathway for generating protein diversity. One noncanonical C-to-U coding edit of interest caused a P124L amino acid change in the activity-regulated Arc1 mRNA that encodes a retroviral-like Gag protein involved in synaptic plasticity (Hantak et al., 2021). Approximately 50% of total Arc1 mRNA was edited at this site in both Ib and Is neurons. Given Drosophila Arc1 assembles into higher order viral-like capsids similar to mammalian Arc proteins (Erlendsson et al., 2020), it will be interesting to determine if and how this amino acid change might alter capsid formation or structure. In terms of the underlying editing machinery, the mammalian APOBEC cytosine deaminases generate C-to-U edits in target mRNAs like apoB (Salter et al., 2016). The Drosophila genome encodes three predicted cytosine deaminases (CG8349, CG8353 and CG8360), though their role in editing has not been studied. Our Patch-seq data indicate CG8353 and CG8360 are expressed in 1b and 1s MNs at similar or higher levels as Adar (Ib TPM (transcripts per million): Adar – 1.51; CG8353 – 1.53; CG8360 –4.77; Is TPM: Adar – 2.16; CG8353 – 1.44; CG8360 – 4.36), suggesting C-to-U editing may be quite robust in these neuronal subtypes. However, given the large number of silent changes to protein sequence, it seems unlikely that noncanonical editing plays a major role in regulating larval MN function compared to ADAR-dependent editing.

Materials and Methods

Single neuron Patch-seq RNA sequencing

RNA sequencing data from 105 MN1-Ib and 101 MNISN-Is MNs (GEO accession #GSE222976) (Jetti et al., 2023) was used for editing analysis. 10XUAS-mCD8-GFP (RRID:BDSC_32185) was expressed using either a Ib-specific (RRID:BDSC_40701) or Is-specific (RRID:BDSC_49227) Gal4 driver. The nucleus and surrounding cytoplasm of GFP+ MNs in 3^rd instar larvae were extracted using a patch pipet and immediately frozen on dry ice. cDNA was generated, purified, dual indexed, and run on an Illumina NextSeq500 (paired-end 75 nt reads) sequencer as described (Jetti et al., 2023).

Genome sequencing

High purity genomic DNA from parental Ib-Gal4, Is-Gal4, and 10XUAS-mCD8-GFP lines was collected and sequenced. Flies were grown at 25°C. 10 adult male flies from each genotype were collected under CO₂ anesthetization into a 1.5mL centrifuge tube and immediately frozen in liquid nitrogen. 120 μL of grinding buffer (0.2 M sucrose, 0.1 M Tris pH 9.2, 50 mM EDTA, 0.5% SDS) was added and a VWR pellet mixer (VWR, cat# 47747-370) was used to homogenize the flies for one minute. Samples were vortexed for 10 seconds and placed on ice for 15 minutes. 500 μL of Phenol:Chloroform:Isoamyl Alcohol (Millipore-Sigma, SKU# 6805-100ML, Omnipur 25:24:1 and TE Buffered Saturated pH 6.7/8.0) was added and vortexed for 30 seconds. Samples were centrifuged at 16,100 relative centrifugal force (RCF) for 5 minutes to pellet debris and separate protein to the interphase and DNA in the upper phase. The DNA-containing upper phase was pipetted into 1 mL of 100% ethanol in a new 1.5 mL tube and vortexed briefly to mix. After incubation overnight at 4°C, samples were centrifuged at 16,100 RCF for 10 minutes to pellet DNA. The pellet was washed in 70% ethanol by pipetting, centrifuged at 16,100 RCF for 10 minutes, and the pellet air-dried for 30 minutes. 50 μL of DNAse-free ultrapure water was added to each sample and incubated at 37°C for 15 minutes before resuspension. DNA samples were prepared as libraries for sequencing using the Nextera DNA Flex (Illumina) kit and sequenced on an Illumina NextSeq500.

RNA edit detection and parental DNA variation calling

Paired-end fastq files were mapped to Drosophila_ BDGP6.92 reference using BWA. Sequencing reads with mapping quality score less than 5 were removed from downstream analyses using Samtools. PCR duplicates were removed using GATK MarkDuplicatesSpark. The variations were first called using Samtools mpileup, and the high-quality SNPs (depth>=8 and quality score >=20) and indels were selected as the database of known polymorphic sites for GATK base quality recalibration. Base quality recalibration was accomplished using GATK BaseRecalibrator and ApplyBQSR. Variations were then called by GATK HaplotypeCaller with VCF outputs. Variations were subset to exonic regions only. Only SNPs were kept for downstream filtering. Since genes localize to both DNA strands while RNA sequencing reads are only aligned to the positive strand, A-to-G editing in minus strand genes appear as T-to-C. To facilitate data analysis, edits in minus strand genes are listed as their reverse complement in the text.

Filtering RNA variations against parental DNA SNPs

SNPs in RNA samples were excluded if they existed in any of the three parental DNA samples as determined by GATK Haplotypecaller and Samtools mpileup based on genomic location and genotype. To prevent mapping complication of nearby indels, RNA SNPs with indels within 10 bp were filtered out. The remaining RNA SNPs were selected as candidate RNA editing sites for further validation.

Filtering of candidate RNA edits

A custom MATLAB code was used to further analyze and filter variant call output from Samtools and GATK. Potential edits were removed if (1) there were ≤ 10 RNA sequencing reads covering this genomic location, (2) the GATK Phred-scaled probability (QUAL score) of a polymorphism was ≤ 20, (3) the cell-specific edit site had ≤ 10% editing rate, (4) the edit was found in fewer than 10 Ib or Is cells, or (5) less than 96% of the sequenced genomic DNA nucleotide at this site failed to match the reference genome (a genomic SNP). Remaining edits were processed to categorize where in the gene the edit occurred. In rare occasions, edits occurred in overlapping genes. In these cases, only edits occurring within exons were kept and listed as an edit in each gene if more than one gene had an exon at this location. Manual inspection revealed many noncanonical edits/SNPs clustered near the beginning of the 5’UTR (where sequencing errors were more common), so edits < 20 nucleotides from the 5’UTR start site were excluded. The remaining high confidence edit sites were analyzed to provide summary statistics, including average editing rate across Ib and Is neurons. Gene boundaries for CDS, 5’UTR and 3’UTR, previously annotated A-to-I editing events, gene descriptions, Pfam locations, and other genome-related information were downloaded from FlyBase on Sept 15 2023.

Software and code

Software versions were as follows: Star (2.5.3a), BWA (0.7.12), Samtools (1.10), GATK (4.1.2.0), MATLAB (2022b). Three-dimensional protein structures were produced using AlphaFold2 (Jumper et al., 2021). The full MATLAB code with instructions for RNA editing analysis is downloadable as a ZIP file from Supplemental Material.

RNA editing conservation analysis

Canonical RNA edits that led to a missense change in the translated amino acid were analyzed to determine if the RNA editing site changed an amino acid that was conserved across species. NCBI BLAST was used to locate the relevant amino acid in polypeptide strings of other insects or mammals and compared to the Drosophila reference genome. Sites were considered conserved if the unedited amino acid at this location was identical between Drosophila and other species.

Comparison of Patch-seq edits with publicly available RNAseq datasets

Raw fastq files were downloaded from existing Drosophila RNA sequencing datasets (GEO Accession ID GSE250462: GSM7979414-7979425 (Cypranowska et al., 2024); GSE37232: GSM914102, 914107, 914108, 914113, 915213-915215, 915220, 914114-914119, 914095 (Rodriguez et al., 2012); GSE113663: GSM3110802-3110804 (Sapiro et al., 2019); GSE223345: GSM6945880-6945885 (Bukhari et al., 2024); EMBL-EBI ArrayExpress E-MTAB-10065: Samples 13-25 EH1_1, EH1_2, EH1_2.2, EH1_3, EH1_3.2, EH1_4, EH1_5, EH1_6, EH1_7, EH1_8, EH1_9, EH1_10, EH2_1, only “Lmerged” files (Harrell et al., 2021) and mapped to Drosophila_melanogaster.BDGP6.32.108 reference genome using STAR (2.7.11b). Neither FastqQC nor adapter trimming steps were performed. The --outSAMtype BAM SortedByCoordinate command was used to output BAM files sorted by coordinate following alignment. The output files “Aligned.sortedByCoord.out.bam” were then indexed by running igvtools within Integrative Genomics Viewer (IGV, version 2.15.2, (Robinson et al., 2011) for visualization. For each genomic dataset loaded into IGV, the genomic coordinate (“edit position”) for each canonical CDS, canonical UTR, and noncanonical CDS editing event from Patchseq was inputted. The number of reads with the edit variant (“ALT”) and those with the unedited variant (“REF”) were recorded. Editing frequency was calculated as ALT/(REF+ALT).

Patch-seq edit positions were also checked in the gDNA_Fly_yw dataset (GSM914095 from GSE37232,(Rodriguez et al., 2012). Edit sites that did not show 100% unedited “REF” reads were eliminated to filter out any additional SNPs. Editing frequencies for the remaining sites were pooled to obtain compiled averages as described below. Editing frequencies recorded for six OK6-GAL4>mCherry-NLS replicates (GSM7979414-7979419, from GSE250462) were averaged to obtain the compiled “Larval” data. Editing frequencies for three OK6-GAL4>mCherry-NLS, RBP-RNAi datasets (GSM7979420-7979422) and three OK6-GAL4>mCherry-NLS, Unc13-RNAi datasets (GSM7979423-7979425) were averaged to obtain RBP and Unc13 knockdown data. Three UAS-unc84-2XGFP; elav-GAL4 datasets (GSM3110802-3110804) were averaged to obtain the compiled “Adult1” data. pAseq_Fly_ZT2_yw (GSM915213), pAseq_Fly_ZT14_yw (GSM915214), pAseq_Fly_ZT0_PE_CS (GM915215), and pAseq_Fly_ZT20_PE_CS (GSM915220) were averaged to obtain the compiled “Adult2” polyA mRNA data. Nascentseq_Fly_ZT2_R1_yw (GSM914102), Nascentseq_Fly_ZT22_R1_yw (GSM914107), Nascentseq_Fly_ZT2_R2_yw (GSM914108), and Nascentseq_Fly_ZT22_R2_yw (GSM914113) were averaged to obtain the compiled “Adult2” Nascent mRNA data. Nascentseq_Fly_FM7a_R1 (GSM914116) and Nascentseq_Fly_FM7a_R2 (GSM914117) were averaged to obtain the compiled “Adult2” control data for comparison to the Adar null mutants (averaged Nascentseq_Fly_ADAR0_R1 (GSM914118) and Nascentseq_Fly_ADAR0_R2 (GSM914119). These compiled averages were used for generating all correlation plots. For all correlation plots, simple linear regression was performed using GraphPad Prism version 10.5.0 (673) for macOS (GraphPad Software, Boston, Massachusetts USA, www.graphpad.com). In all cases, regression lines were not constrained; the lines were not forced to pass through the origin (0,0).

For supplemental tables 8 and 9, each site with a compiled average of at least 10% editing in one of the three datasets was listed as confirmed. Any site with less than 10% editing was listed as not confirmed. Edits with no read coverage were listed as no reads. The number of canonical CDS, canonical UTR, and noncanonical CDS edits belonging to each dataset (Larva, Adult1 or Adult2) were determined and plotted to generate pie charts.

Statistical analysis

For comparisons between two groups, statistical significance was determined using an unpaired, two-tailed, two-sample Student’s t-test. For comparisons between three or more groups, statistical significance was determined using a one-way ANOVA followed by the Tukey-Kramer post-hoc multiple comparisons test. Figures depict the mean of each distribution and individual data points, and asterisks denote p values of: *, p≤0.05; **, p≤0.01; ***, p≤0.001; and ****, p≤0.0001. Normalcy of data distribution was determined using the one-sample Kolmogorov-Smirnov test. The correlation coefficient (r) between two sets of data was determined using the Pearson correlation method.

Figure supplements

Figure 7 – figure supplement 1.

Figure 7 – figure supplement 2.

Figure 8 – figure supplement 1.

Figure 8 – figure supplement 2.

Figure 8 – figure supplement 3.

Acknowledgements

This work was supported by NIH grants MH104536 and NS117588 to J.T.L. AB.C. was supported in part by NIH pre-doctoral training grant F31NS118948. We thank the Bloomington Drosophila Stock Center (NIH P40OD018537) for providing Drosophila strains, Duanduan Ma and Charlie Whittaker for help with genomic and editing analysis, and members of the Littleton lab for helpful discussions and comments on the project and manuscript.

Additional information

Data Availability

The source data supporting this study are available within the paper, supplementary material and indicated repositories. RNA sequencing data is available at the GEO data repository (accession #GSE222976). DNA sequencing data for parental genotypes is available at the NCBI BioSample Database (accession #PRJNA1112347). MATLAB codes are provided in Supplemental Material.

Author contributions

Conceptualization, A.B.C and J.T.L.; methodology, A.B.C., M.O.I. and S.K.J.; investigation, A.B.C., M.O.I. and J.T.L.; funding acquisition, J.T.L.; supervision, J.T.L.; writing—original draft, A.B.C., M.O.I. and J.T.L.; writing—review & editing, A.B.C., M.O.I., S.K.J. and J.T.L.

Funding

National Institute of Mental Health (MH104536)

National Institute of Neurological Disorders and Stroke (NS117588)

Additional files

Supplemental Table 1. Canonical RNA (A-to-I) editing sites identified across all individual Ib MNs, Is MNs and abdominal muscles of 3^rd instar larvae.

Supplemental Table 2. Summary of canonical RNA editing sites identified in Ib MNs, Is MNs and abdominal muscles.

Supplemental Table 3. Noncanonical RNA (non-A-to-I) editing sites identified across all individual Ib MNs, Is MNs and abdominal muscles of 3^rd instar larvae.

Supplemental Table 4. Summary of noncanonical RNA editing sites identified in Ib MNs, Is MNs and abdominal muscles.

Supplemental Table 5. Canonical RNA editing targets with the highest editing rate that cause amino acid substitutions.

Supplemental Table 6. Canonical RNA editing site differences in Ib versus Is MNs following removal of differentially expressed genes (DEGs).

Supplemental Table 7. Canonical RNA editing sites identified in larval abdominal muscles.

Supplemental Table 8. Comparison of canonical RNA editing sites identified via Patch-seq to other RNAseq datasets.

Supplemental Table 9. Comparison of noncanonical RNA editing sites identified via Patch-seq to other RNAseq datasets.

MATLAB code. Zip file of MATLAB code for RNA editing detection. This file contains the necessary MATLAB codes to process output of GATK/SAMtools variant callers as described in methods.