Characterization of high-confidence A-to-I editing sites across retina development

To obtain a global landscape of A-to-I editing in the developing mouse retina, we performed total RNA-seq to profile five time points, including E13, P0, P6, P21, and P42. In total, the RNA-seq experiments yielded more than 1.19 billion raw reads, with an average of ~119 million reads per sample (Fig. 1A; Additional file 1: Supplementary Table 1). After quality control and data preprocessing, we applied REDItools2 [24] to characterize the RNA editing profiles. Our subsequent filtering steps, illustrated in Fig. 1A and Additional file 2: Fig. S1, resulted in a set of 17,874 high-confidence RNA editing sites, of which 15,109 were A-to-I editing sites, making up 84.53% of the entire set (Fig. 1B; Additional file 1: Supplementary Tables 2 and 3). This proportion is in line with prior research findings [25, 26]. Among these A-to-I editing sites, 8104 were previously reported by REDIportal, while the remaining 7005 were newly discovered, and both exhibited a similar motif pattern (Additional file 2: Fig. S2 and Fig. S3). Sanger sequencing of cDNA and gDNA further validated some newly discovered sites, including two sites specific to the mouse retina within Rgs9bp, thereby confirming their authenticity (see the “Methods” section; Additional file 2: Fig. S4). Moreover, the newly discovered sites were found to have a significant association with functions such as synapse and visual perception (Additional file 2: Fig. S3). At least, these results indicated that some, if not all, of these newly discovered sites are retina-specific. Given that A-to-I editing is the most prevalent type of editing and has significant impacts on development [14, 27, 28], we chose to focus our subsequent analysis solely on this type of editing.

Fig. 1

Genome-wide characterization of high-confidence A-to-I editing sites. A Schematic illustration of the experimental design and high-confidence A-to-I editing site identification and annotation. B Distribution of RNA editing types. The bar graph displays the number of each type of RNA editing, with “AG” representing A-to-I editing sites. C Nucleotide context around A-to-I editing sites, consistent with previous findings. D Heatmap of editing levels in different developmental time points. A value closer to 1 indicates a higher similarity of editing levels between samples. E Principal component analysis of editing levels in different developmental time points. F Distribution of editing sites on genomic regions. The inserted pie displays their relative fraction. Nonsynonymous refers to editing sites in the CDS that result in changes in amino acids, while synonymous refers to editing sites that do not cause such changes

We found that the level of guanosine was lower in the nucleotide before the editing site and higher in the nucleotide after, which aligns with the substrate requirements of ADAR editing (Fig. 1C). Hierarchical clustering analysis showed high consistency in editing levels between replicates and clear separation between developmental time points (Fig. 1D). The results of the principal component analysis reflected a developmental progression, from the embryonic (E13) to neonatal (P0 and P6) and then to eye-opening (P21 and P42), in line with the maturation of the retina (Fig. 1E). Consistent with prior studies [28,29,30], we found that the majority (59.9%) of A-to-I editing sites were presented within introns, with only a small proportion (0.9%) located within coding sequences (CDS). Among those within CDS regions, 73.53% led to non-synonymous changes (Fig. 1F). Overall, our analysis demonstrates the high reliability of the A-to-I editing sites we identified.

Temporal changes of A-to-I editing across retina development

To investigate temporal dynamics of A-to-I editing, we initially quantified editing sites and noted substantial variation across time points, ranging from 310 to 11,014 sites (Fig. 2A). Despite this variability, a marked increase in site numbers was observed over time, with a steep surge occurring post eye-opening. Notably, this phenomenon could not be attributed to sequencing coverage bias, as no significant correlation between the number of A-to-I editing sites and sequencing depth was observed (Additional file 2: Fig. S5). We further categorized these editing sites into five groups based on their developmental prevalence. The majority (67.38%, 10,181 sites) were exclusive to a single time point, with only a negligible fraction (0.66%, 99 sites) shared across all time points (Fig. 2B). Analysis between different groups revealed substantial disparities in editing levels, with a trend of higher editing levels as prevalence increased (Fig. 2C). Intriguingly, shared editing sites displayed elevated editing activity that increased during retina development (Fig. 2C). These sites preferentially localized to 3′ UTR regions (Fig. 2D; Additional file 2: Fig. S6). Our enrichment analysis revealed that they were often located within genes associated with functions such as “regulation of mRNA processing,” “ATP-dependent chromatin remodeling,” and “RNA splicing” (Fig. 2E), suggesting a possible functional purpose for these sites. In contrast, timepoint-specific editing sites displayed relatively low levels of editing activity, which might be due to purifying selection that impedes their editing ability or prevalence [31].

Fig. 2

Temporal changes of A-to-I editing. A Temporal distribution of A-to-I editing sites. B Intersection of editing sites at each time point. C Editing levels grouped by prevalence. The “f1” refers to editing sites that appeared only once in the 5 time points, while “f5” refers to editing sites that are shared by all 5 time points. The editing levels of f5 in different developmental time points were shown on the right. Significance testing was performed by the Wilcoxon rank sum test, and p-values were shown at the top of the box plots. D Proportion of sites in different genomic regions concerning prevalence. It shows the proportion of sites in 3′-UTR increasing and the proportion of sites in intronic decreasing with increasing prevalence. E Enrichment analysis of editing genes shared by all developmental time points. The top-ranked enriched GO terms are shown here. F Bar plot showing the number of differential editing sites between adjacent developmental time points. G Enriched GO biological process terms of genes with differential up- and downregulated editing sites. H Normalized expression of three ADARs at the transcriptional (left) and translational level (right). I Pearson’s correlation analysis between the editing number of all sites and ADAR expression at the transcriptional and translational levels. J Pearson’s correlation analysis between editing levels of completely shared sites and ADAR expression at the transcriptional and translational levels, respectively

The temporal changes of A-to-I editing were then studied by analyzing the differences in editing levels between adjacent time points using REDITs [32]. Our results showed that a total of 604 sites underwent differential editing, with a remarkable transformation in the number of differentially edited sites before and after eye-opening (Fig. 2F and Additional file 2: Fig. S7; Additional file 1: Supplementary Table 4). Our enrichment analysis indicated that as the retina developed, editing sites that experienced upregulation were primarily linked to synaptic vesicles, such as “vesicle−mediated transport in synapse” and “synaptic vesicle transport,” while those that were downregulated were mainly related to RNA splicing, such as “mRNA splicing, via spliceosome” and “regulation of mRNA splicing, via spliceosome” (Fig. 2G). These results suggest the involvement of A-to-I editing in retina development, particularly in modulating neurotransmitter release from synaptic vesicles and guiding alternative splicing decisions.

Furthermore, we investigated the relationship between ADAR expression and RNA editing. We found that ADAR2 had the most noticeable increase in expression throughout retina development, as seen at both transcriptional and translational levels, compared to the other two ADAR genes (Fig. 2H). Additionally, we found a significant positive relationship between the number of editing sites and ADAR2 expression, as indicated by Pearson’s correlation coefficients of 0.92 (p-value = 0.027) and 0.98 (p-value = 0.003) for transcription and translation, respectively (Fig. 2I). While no clear correlation was found between editing levels of all editing sites and ADAR1/2 expression, a marginally significant negative correlation emerged with the translational level of ADAR3 (Additional file 2: Fig. S8). Despite this, a portion of the editing variability could be attributed to ADAR expression, as indicated by a significant positive relationship between its transcription and editing levels of editing sites completely overlapped across all time points (Pearson’s r = 0.91, p-value = 0.034) and a marginally significant positive relationship between its translation and editing levels at the same sites (Pearson’s r = 0.87, p-value = 0.054) (Fig. 2J). Specifically, we found that the identified A-to-I editing sites covered approximately 46% of previously reported ADAR2 targets [12], but only 10% of ADAR1 targets (Additional file 2: Fig. S9). Collectively, our results suggest that ADAR2 may play a more important role without exclusive regulation by ADAR1/3 on RNA editing.

Timepoint-specific A-to-I editing pattern on retina development

We next explored the RNA editome in greater detail to understand the changes in the editing pattern. By using mfuzz clustering [33], we identified six distinct groups of temporal editing profiles, as shown in Fig. 3A and Additional file 1: Supplementary Tables 5 and 6. The first cluster (c1), consisting of 1885 editing sites, showed a pattern of concurrent editing, with a sudden increase in editing levels following eye-opening. This pattern was also observed in cluster 2 (c2), which was made up of 1956 editing sites (Fig. 3B). Further analysis revealed that the sites within c1 and c2 were associated with functions such as “regulation of long-term neuronal synaptic plasticity” and “sensory perception of light stimulus” (Fig. 3C; Additional file 1: Supplementary Table 7). The editing patterns of clusters 3–6 were unique to their respective time points. Cluster 3 (c3), which was comprised of editing sites specific to P0, was characterized by functions related to those such as “regulation of DNA metabolic process,” “DNA repair,” and “covalent chromatin modification,” coinciding with actively cellular differentiation at this time point, marked by the formation of a substantial number of rod cells [34]. Interestingly, Crx, a crucial transcription for photoreceptor cell differentiation, was also edited during this time point. Cluster 4 (c4), made up of P6-specific editing sites, was characterized by functions related to those such as “neuron projection arborization” and “negative regulation of binding.” Cluster 5 (c5), made up of P21-specific editing sites, was characterized by functions related to those such as “synapse organization” and “vesicle−mediated transport in synapse,” and cluster 6 (c6), made up of P42-specific editing sites, was characterized by functions related to such as “covalent chromatin modification” and “histone modification,” suggesting that proper editing of sites in the retina may be necessary for functional maturation of retinal cells.

Fig. 3

Developmental patterns of A-to-I editing. A Heatmap displaying the temporal distribution of A-to-I editing during development, with the right showing the distribution of retinal markers in these patterns. The line plot depicts the editing level trends for each pattern, with colored lines indicating the normalized mean editing level of each pattern, and gray lines representing the normalized editing level of an individual editing site. B Bar plots showing the number of genes and editing sites included in each editing pattern. C The top-ranked 5 or all enriched GO terms for each pattern in panel A are listed, with retina markers shown in parentheses after each function category

Interplay between alternative splicing and A-to-I editing

To gain insights into the connection between RNA editing and alternative splicing, we examined their interrelation. By analyzing transcriptome data (see the “Methods” section; Additional file 1: Supplementary Table 8), we found that 73% of genes with RNA editing also exhibited alternative splicing (Fig. 4A). The presence of RNA editing was found to be significantly more prevalent among genes with splicing compared to those without, regardless of whether the transcript length was normalized (Fig. 4B; Additional file 2: Fig. S10; Fisher’s exact test, p-value < 0.01). When comparing genes with and without RNA editing but both with alternative splicing, we observed that edited genes had a greater number of alternative splicing events per gene (Fig. 4C) and higher splicing efficiency, quantified by the percent spliced in (PSI) value (Wilcoxon rank sum test, p-value < 2.2e−16) (Fig. 4D). Our analysis also revealed a proximity between RNA editing and intron-retained (IR) events (Additional file 2: Fig. S11), suggesting that RNA editing may have a greater impact on IR events compared to other events, such as exon skipping and mutually exclusive exons (EX), alternative acceptors (Alt3), alternative donors (Alt5), and exon skipping for micro exons (MIC).

Fig. 4

Interplay between A-to-I editing and alternative splicing. A Proportion and number of editing genes with and without alternative splicing in each developmental time point. B Association between A-to-I RNA editing and alternative splicing, by classifying genes into four categories based on the presence or absence of editing sites and alternative splicing events, and using Fisher’s exact test to determine significance (p-value < 0.01). C Comparison of the average number of splicing events between genes with and without editing. D Comparison of the percent spliced in (PSI) values between editing and non-editing genes. E Patterns of Normalized PSI values (left) and normalized editing levels (right) for strongly correlated pairs of splicing events and nearby editing sites. The colored line represents the overall trend of changes for each group, while each gray line reflects the change of a specific splicing or editing event. F Heatmap displaying the tendency of editing levels and splicing efficiency for strongly correlated pairs of splicing events and editing sites, with “pos” and “neg” indicating positive and negative correlation, respectively. G Heatmap displaying the top-ranked 5 or all enriched GO terms for six strongly correlated groups. H Comparison of the distances between paired positively correlated splicing events and editing sites with those of negatively correlated pairs. I Mosaic plot showing the distribution of different types of splicing events between positively and negatively correlated pairs of splicing events and editing sites. Significance testing was performed by Wilcoxon rank sum test: ns, not significant; ****p < 0.0001

To further examine the extent of their developmental interrelation, we focused on analyzing 9390 pairs of editing sites (excluding those near the 4 nt intronic side of the splicing sites) and their corresponding nearby splicing events. We found that 2143 pairs were strongly correlated and categorized them into six distinct groups using Mfuzz (see the “Methods” section; Fig. 4E; Additional file 1: Supplementary Table 9). Our results showed that the changes in editing level and splicing efficiency followed a similar trajectory in groups 2, 3, 5, and 6, while groups 1 and 4 displayed a contrasting trajectory. Enrichment analysis showed that functions related to chromosome and mitochondrion, such as “positive regulation of chromosome organization” and “mitochondrion disassembly,” were over-represented in groups 5 and 6. This suggests that editing and splicing were closely intertwined in shaping message RNA. On the other hand, functions related to development and modification, such as “dendrite development” and “tRNA modification,” were over-represented in groups 1 and 4 (Fig. 4G). Correlation analysis between editing level and splicing efficiency showed that 957 pairs of editing sites and splicing events had a significantly strong relationship (absolute Pearson’s r ≥ 0.7 and p-value ≤ 0.05) (Fig. 4F). Positively correlated editing sites and splicing events were located closer together than negatively correlated ones (Wilcoxon rank sum test, p-value = 0.00066) (Fig. 4H and Fig. S7). IR events were more frequent than expected by chance in positive correlations, while EX events were dominant in negative correlations, suggesting that the impact of RNA editing may vary depending on the type of splicing event, with high editing activity tending to favor the preservation of nearby intron or the suppression of nearby exon (Fig. 4I).

Alteration of translatome conferred by A-to-I editing

Translation rate and output can be impacted by RNA editing, but the extent of this impact is yet to be determined. To shed light on this, we used ribosome profiling to generate translation profiles, yielding an average of ~66 million raw reads per sample (see the “Methods” section; Additional file 1: Supplementary Table 1). Our results revealed that the combination of alternative splicing and RNA editing (AS & ES) resulted in the highest average number of actively translated transcripts per gene (Fig. 5A). This was followed by the group with only splicing (AS & Non-ES), then by the group with only editing (Non-AS & ES), and finally by the group with neither editing nor splicing (Non-AS & Non-ES). However, the Non-AS & Non-ES group had the highest translational efficiency, followed by the Non−AS & ES, AS & Non-ES, and AS & ES groups (Fig. 5B). These findings indicate that editing and splicing can increase coding capacity and diversify the translatome, with a synergistic effect when used together, for example, the editing level and splicing efficiency of retina-specific gene Pcdh15, Stx3, Pde6b, and Tia1 synergistically inhibit their TE (Additional file 2: Fig. S12). Notably, the increasing diversity of translated transcripts was accompanied by heightened usage of RNA editing and splicing, with the latter having a more pronounced impact on gene’s translational efficiency than the former.

Fig. 5

A-to-I editing induces changes in translatome. A Point-range plot displaying the mean number of translatable transcripts for four groups of genes (AS & ES, genes with both alternative splicing events and RNA editing sites; AS & Non-ES, genes with alternative splicing events but none editing sites; Non-AS & ES, genes without alternative splicing events but have editing sites; Non-AS & Non-ES, genes without either alternative splicing events or editing sites) classified based on the presence or absence of splicing events or editing sites. B Boxplot showing translational efficiency between different groups of genes (the groups classified as same as in A). C Number of genes with differential splicing events, differential editing sites, and differential translational efficiency in pairwise comparisons. D Violin plot comparing translational efficiency among four gene groups classified based on whether they exhibited differential editing or differential splicing efficiency. E Heatmap showing the top-ranked 5 enriched biological process GO terms for different gene groups in panel D. Significance testing was performed by Wilcoxon rank sum test: ****p < 0.0001

In light of these findings, we further examined the effect of RNA editing on gene translation. We performed differential translational efficiency (dTE) analysis between adjacent time points and found 2936 dTE genes that mirrored retina development, with two pronounced peaks in gene number between E13 and P0, and P6 and P21 (see the “Methods” section and Fig. 5C; Additional file 1: Supplementary Table 10). In parallel, we also found 3453 differential splicing efficiency (dPSI) genes and 440 differential editing level (dEL) genes (see the “Methods” section). When dTE genes were classified into four groups based on their dPSI or dEL status, we observed that in the Non-dAS & Non-dES group, there existed a close balance between up- and downregulated dTE genes, with 51.50% and 48.50%, respectively. The balance was disrupted in the presence of dPSI or dEL, resulting in the majority of dTE genes being downregulated in the Non-dAS & dES (65.98%), dAS & Non-dES (67.91%), and dAS & dES (88.1%) groups (Fig. 5D). These results suggest that RNA editing and splicing serve as a buffering mechanism to reduce gene translational efficiency, with both having a coordinated effect. Enrichment analysis further revealed that only the Non-dAS & Non-dES group had an over-representation of functions related to retina development, such as “visual perception” and “visual system development,” while the other three groups had an over-representation of functions related to the basic processes of life, such as “chromatin silencing” and “regulation of chromosome organization” (Fig. 5E).