
==== Front
RNA Biol
RNA Biol
RNA Biology
1547-6286
1555-8584
Taylor & Francis

39256954
10.1080/15476286.2024.2397757
2397757
Version of Record
Research Article
Research Paper
An orthology-based methodology as a complementary approach to retrieve evolutionarily conserved A-to-I RNA editing sites
J. LIU ET AL.
RNA BIOLOGY
Liu Jiyao #
Zhao Tianyou #
Zheng Caiqing #
Ma Ling
Song Fan
Tian Li
Cai Wanzhi
Li Hu
https://orcid.org/0000-0003-2311-9859
Duan Yuange
Department of Entomology and MOA Key Lab of Pest Monitoring and Green Management, College of Plant Protection, China Agricultural University , Beijing, China
CONTACT Yuange Duan duanyuange@cau.edu.cn Department of Entomology and MOA Key Lab of Pest Monitoring and Green Management, College of Plant Protection, China Agricultural University, Beijing 100193, China
# Co-first author

10 9 2024
2024
10 9 2024
21 1 2945
Integra10 9 2024
Integra10 9 2024
14 8 2024
19 8 2024
© 2024 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
2024
The Author(s)
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.

ABSTRACT

Adar-mediated adenosine-to-inosine (A-to-I) mRNA editing is a conserved mechanism that exerts diverse regulatory functions during the development, evolution, and adaptation of metazoans. The accurate detection of RNA editing sites helps us understand their biological significance. In this work, with an improved genome assembly of honeybee (Apis mellifera), we used a new orthology-based methodology to complement the traditional pipeline of (de novo) RNA editing detection. Compared to the outcome of traditional pipeline, we retrieved many novel editing sites in CDS that are deeply conserved between honeybee and other distantly related insects. The newly retrieved sites were missed by the traditional de novo identification due to the stringent criteria for controlling false-positive rate. Caste-specific editing sites are identified, including an Ile>Met auto-recoding site in Adar. This recoding was even conserved between honeybee and bumblebee, suggesting its putative regulatory role in shaping the phenotypic plasticity of eusocial Hymenoptera. In summary, we proposed a complementary approach to the traditional pipeline and retrieved several previously unnoticed CDS editing sites. From both technical and biological aspects, our works facilitate future researches on finding the functional editing sites and advance our understanding on the connection between RNA editing and the great phenotypic diversity of organisms.

KEYWORDS

A-to-I RNA editing
honeybee
new methodology
novel editing sites
caste-specific editing
National Natural Science Foundation of China 10.13039/501100001809 32300371 the Young Elite Scientist Sponsorship Program by CAST the Young Elite Scientist Sponsorship Program by BAST the 2115 Talent Development Program of China Agricultural University This study is financially supported by the National Natural Science Foundation of China [no. 32300371] the Young Elite Scientist Sponsorship Program by CAST [no. 2023QNRC001], the Young Elite Scientist Sponsorship Program by BAST [no. BYESS2023160], and the 2115 Talent Development Program of China Agricultural University.
==== Body
pmcIntroduction

A-to-I RNA editing in organisms

The ADAR (adenosine deaminase acting on RNA) protein family mediates the adenosine-to-inosine (A-to-I) mRNA editing, which is the most abundant and prevalent RNA modification in metazoans [1–3], fungi [4,5], and bacteria [6,7]. In different species, thousands to millions of adenosine sites in the mRNAs are potentially editable [8,9]. Due to the structural similarity between inosine and guanosine, A-to-I RNA editing causes similar effects to A-to-G DNA mutations [10–12], potentially changing the genomically encoded amino acids (AA), leading to a phenomenon called ‘recoding’ which simply means nonsynonymous RNA editing [13,14]. An intuitive advantage of RNA recoding comes to its flexibility to adjust the proteomic diversity, the theory of which was summarized as the ‘proteomic diversifying hypothesis’ [15]. This adaptive hypothesis has been experimentally verified in a few species like Fusarium graminearum (fungi) [16] and cephalopods (mollusc) [17,18]. For example, for a conserved recoding site in Fusarium graminearum, the pre-editing isoform is fitter under asexual stage and the post-edited isoform is fitter under sexual stage, and RNA editing can choose which allele to be used under a specific condition, increasing the overall fitness [16]. In cephalopods, the post-edited motor protein has higher motility than the pre-editing isoform under cold temperature, serving as a compensatory mechanism for the animal [17,18]. Another highly conserved recoding site in potassium channel also have similar effects [19,20]. These fascinating cases all highlight the advantage of having the editing ability compared to a hardwired DNA allele. The continuous finding of adaptive RNA editing sites will advance our understanding of the environmental adaptation, molecular and phenotypic plasticity, and epigenetic regulation of organisms.

However, to perform functional validation, it is quite challenging to select such adaptive recoding sites from the sea of total RNA editing events. One potential approach to narrow down the candidate functional sites is to look at the conserved recoding sites across multiple species [21,22]. But it is important to note that the prerequisite for finding conserved editing sites is to first identify the reliable lists of editing sites in each target species, respectively. Here, we summarize the widely recognized challenges in the identification of RNA editing sites, and then proposed a novel complementary methodology to retrieve potential conserved recoding sites which might be missed by the traditional RNA editing detection pipeline.

Challenges in the identification of RNA editing sites

In most species with available reference genome and transcriptome, the accurate identification of a reliable list of RNA editing sites is still challenging [23]. First, one needs to distinguish bona fide RNA editing sites from the sequencing errors in the transcriptome data. The sequencing error rate of next generation sequencing (NGS) is around 0.1–1%, and only those RNA editing sites with editing levels much higher than 1% can be reliably identified. However, most recoding sites have editing level lower than 10% in nervous systems [20] and the use of other tissues or whole body (like small insects) will further dilute the editing signal [23]. Second, another major confounding factor is the single nucleotide polymorphism (SNP) in the genome. The mismatch (MM) between RNA-Seq and the reference genome is not necessarily RNA editing events. Thus, usually a matched DNA-resequencing library from the same individual(s) is needed [23]. But the genomic regions with insufficient DNA-Seq coverage will decrease the confidence of RNA editing detection, let alone in many cases (like non-model insects), a matched DNA-Seq from the same individual is not feasible at all. Third, even if the difficulties in sequencing issues are resolved, the downstream bioinformatic analysis is also full of pitfalls. The mis-alignments between sequencing reads and the reference genome will cause numerous false positive mismatches that mislead the identification of RNA editing sites [23]. Mis-alignments are inevitable due to the existence of repetitive elements, paralogs, and splicing junctions [24]. A way to partially reduce mis-alignments is to require a read to be mapped to a single genomic location [25], but this stringency might affect the reads distribution and distort the genuine editing levels. Other commonly used stringent criteria include high cut-offs on sequencing coverage, high cut-offs on editing level, and requiring a variation site to appear in multiple RNA-Seq samples. These stringent criteria might preclude some true-positive editing sites: e.g. differential editing sites between samples, or a site with randomly uneven coverages between different samples.

Trade-off between quantity and quality exists for the de novo identification of RNA editing due to the lack of prior knowledge

Given the above-mentioned challenges and limitations, the de novo identification of RNA editing sites is highly sensitive to false-positive sites. Here, ‘de novo identification’ means that we do not have any prior knowledges on which sites are more likely to be RNA editing sites so that every adenosine in the genome should be treated with equally stringent criteria. De novo identification is usually performed in the first RNA editing study of a particular species. Then, ‘false-positive sites’ refers to the scenario where we think a site is A-to-I RNA editing but it is actually a SNP or sequencing error. To reduce false-positive rate, a conceivable effort is to obtain a high A>G% fraction among the total RNA variation sites after a series of filters. For example, in the cephalopod [20] and ant [25] studies, over 90% of the final RNA-DNA differences (RDD) are A-to-G variations, representing reliable A-to-I RNA editing. In contrast, some RNA editing studies with low A>G% fraction might be suspected for high false-positive rates. For example, a study of human RNA editing [26] was commented by other groups [24,27], and a study on viral RNA editing [28] was commented by following papers [29,30]. As we will show in the results, more stringent criteria can reduce the false-positive rate (and increase the A>G%), but the total number of A-to-G sites will also dramatically decrease. This means that many true-positive A-to-I RNA editing sites, especially the lowly edited ones, might be eliminated by the strict criteria. In other words, quality comes at the cost of quantity. This quality-quantity trade-off is inevitable during the de novo identification of RNA editing sites.

However, if we have a prior knowledge of RNA editing sites in other species (e.g. Drosophila), then the identification of conserved RNA editing in this species (e.g. Apis mellifera) can be less stringent. For example, if an A-to-G variation is observed in the transcriptome of honeybee and this position is already a known A-to-I RNA editing site in fruitfly, then it might require less stringent cut-offs to convince that this is a conserved RNA editing site between bee and fly. When sequencing error is excluded by reasonable cut-offs, there is almost no worry about SNPs on this site because there are currently no reports on an A-to-G SNP in one species that corresponds to the known A-to-I RNA editing at the orthologous site of another species; the later scenario is virtually rarer than balancing selection on both species. This merit is useful when a matched DNA-resequencing is not available for an RNA-Seq sample. More importantly, this orthology-based methodology focuses on the already known editing sites in e.g. fly and directly checks the status of their orthologous positions in e.g. bee, then there is not the step of A>G% enrichment. As we will show, this step in the de novo identification usually discards many true-positive sites, and thus the orthology-based methodology might be able to retrieve those sites.

It is essential to note that this orthology-based methodology is only a complementary approach to the de novo identification. It does not mean that the new method alone can identify the full set of RNA editing sites, neither can the orthology-based methodology identify more editing sites than the traditional pipeline. The de novo identification is still the most prevalent way to obtain most of the true-positive editing sites in a species, and the orthology-based method benefits us by adding several previously missed conserved RNA editing sites, making the candidate list more complete.

Eusocial insects, transcriptomic plasticity, and RNA editing

Compared to the broad animal species, eusocial insects have evolved the ability to generate phenotypically differentiated individuals under the same set of genome [31,32]. This raises a possibility that the epigenetic and post-transcriptional regulations might play an essential role in and shaping their diverse phenotype. Among several popular ‘beyond-genome’ regulatory mechanisms, methylation generally controls the expression level of genes while A-to-I RNA editing is able to change the protein sequence. This prompts the demand for understanding (1) whether and how A-to-I RNA recoding contributes to the great plasticity of eusocial insects, and (2) can we identify any highly conserved or adaptive RNA recoding sites in their transcriptome.

Despite the major challenges in de novo RNA editing detection in insects, the genome-wide A-to-I RNA editomes of three eusocial insects have been systematically characterized in the past decade, including a study on leaf-cutting ant Acromyrmex echinatior in 2014 [25], bumblebee Bombus terrestris in 2019 [33], and honeybee Apis mellifera in 2021 by ourselves [34]. Caste-specific RNA editing sites were found in bumblebee and leaf-cutting ant, indicating a potential contribution of RNA editing to caste-differentiation.

However, for several reasons, an older version genome assembly of honeybee (A.mel 4.5, http://hymenopteragenome.org/beebase/) was not optimal for identifying RNA editing sites. This fact is reflected by the low abundance of CDS editing sites identified in our previous study, together with a message from our personal contact with one author of a methodological paper [35]. This limitation prevented us from obtaining an accurate landscape of the conservation, adaptation, and regulation of RNA editing in honeybees. Particularly, due to the omission of some true-positive editing sites in honeybee, a conserved editing site between two species (e.g. between honeybee and bumblebee, or between honeybee and fruitfly) will be mis-regarded as a species-specific editing site.

Aims and scopes

Given the imperfectness of the previously de novo identified A-to-I RNA editome in honeybee Apis mellifera [34], together with the known fact of caste-specific RNA editing in bumblebee [33] and leaf-cutting ant [25], the following questions are eager to be answered. (i) How many potentially conserved RNA editing sites in CDS can be retrieved by the orthology-based methodology? (ii) Can we find conserved editing sites that show caste-specificity in multiple eusocial insects? In this work, using the new genome of honeybee and the orthology-based methodology for identification of RNA editing sites, we provide interesting findings and implications on how the transcriptomic plasticity at molecular level might contribute to the remarkable phenotypic divergence in eusocial insects. The methodology is also helpful in finding conserved editing sites between species, narrowing down the candidate functional sites for future detailed studies.

Results

Repeating the traditional RNA editing pipeline on the new genome version of honeybee

Recently, a new genome version HAv3.1 of Apis mellifera was updated. This genome version has remarkably higher completeness and N50/N90, and lower numbers of scaffolds (Table 1), enabling us to more accurately map the sequencing data to the genome and get a more comprehensive profile of RNA editing sites.Table 1. The statistics of two versions of honeybee genomes.

Version	Genome size (Mb)	GC content	Complete-ness	N50 (Mb)	N90 (Mb)	# of scaffolds	
A. mel 4.5	234.1	32.7%	97.6%	0.99	0.11	5644	
HAv3.1	225.3	32.5%	99.4%	13.6	1.07	177	

To understand how exactly the traditional RNA editing detection pipeline might miss some true-positive RNA editing sites, and to what extent can the orthology-based methodology retrieve the missing sites, we first need to schematically redo a similar data analysis to our previous study [34] on the new honeybee genome. The following processes are conceptual and aim at providing a basic understanding of why the de novo identification of RNA editing sites is over-stringent and can miss a plenty of true-positive sites. The head, thorax, and abdomen transcriptomes from two male drones were used (denoted as individual#1 and individual#2), and the DNA-Seq from the same individuals was also incorporated into the analysis (Materials and Methods). Since A-to-I RNA editing is most abundant in insect heads [34], the head transcriptome was the first to be inspected. Take individual#1 for instance, by mapping the head RNA-Seq to the reference genome, 211,448 genic variations sites were obtained and only 34,011 (16.08%) were A-to-G variations. This A>G% is not ideal and should be improved. By requiring FDR < 0.05 in the binomial test to exclude sequencing errors (Materials and Methods), the A>G% became 16.84% (19065/113221) (Supplementary Figure S1). This A>G% is still close to the random expectation and cannot be regarded as RNA editing signal at all [24,26–30].

One may envision that the variations in RNA-Seq against the reference genome can also contain many SNPs. Then, we required the variation sites to be supported by ≥10 reads in DNA-Seq and no alternative alleles appear in DNA. This procedure ensures that the RNA-Seq variants are authentic RNA-DNA differences in that individual. Expectedly, this step remarkably elevated the A>G% fraction but the number of sites also inevitably decreases (811/1463 = 55.43%, Supplementary Figure S1). To further optimize the A>G% fraction, we further asked a variation to appear in heads of the two individuals, in which both samples meet the above criteria. This requirement again enriched the A-to-G variations (321/417 = 76.98%, Supplementary Figure S1). Since the A>G% fraction is highly acceptable, these 321 A-to-G variations are regarded as A-to-I RNA editing sites in drone heads. Then, thorax and abdomen were analysed with identical procedures, and the identified editing sites were combined with the head sites to obtain a union of 407 sites, representing the A-to-I RNA editome in honeybee produced from the traditional identification pipeline (Supplementary Figure S1).

Among the 407 editing sites from de novo identification, only five sites were conserved between honeybee and D. melanogaster (two recoding and one synonymous site in gene Shab and two recoding sites in gene qvr). Given the prevalence of RNA recoding on neuronal genes in insects [36], one should expect more conserved editing sites to be found. This raises the need to develop a new orthology-based method to retrieve the conserved sites which are discarded by the strict criteria of de novo identification. The new methodology should be understood as a complementary approach rather than a better/independent approach. In the following sections, we will show its validity and efficiency.

Combination of known CDS RNA editing sites in representative species

To perform the orthology-based methodology and see whether this method can complement the traditional pipeline, we retrieved the transcriptomes of multiple tissues and castes of Apis mellifera, including the head/thorax/abdomen transcriptomes of a male drone, and the brain transcriptomes of female workers (six nurses and six foragers) (see Materials and Methods for the data accession IDs and the rationale for sample selection). These samples were used for the identification of A-to-I RNA editing sites (Figure 1A). Our general notion is to refer to the known CDS editing sites in several representative species and see whether editing signal can be detected at the orthologous positions in honeybee transcriptome (Figures 1B,C). Notably, the original drone data contain two individuals with polyA+ RNA-Seq, but the individual#2 is of low quality and thus we only used individual#1 (Materials and Methods). We also downloaded the bumblebee transcriptomes [33] to check if the conserved editing sites between honeybee and fruitfly were also detected in bumblebee. This serves as a supporting evidence to show that the editing sites found in honeybee were not artefacts (Materials and Methods). Figure 1. Data collection and scheme for identifying RNA editing sites. (A) Different transcriptome samples of honeybee used in this study. (B) The blastp procedure to determine the orthologous gene (CDS) pairs between D. melanogaster and honeybee (Materials and Methods). The first round of blastp required reciprocal best and the second round of blastp was done for the remaining genes and required E-value <1E–6 and matching length >50% on both sides. D. melanogaster was shown as an example and the matching between C. chinensis and honeybee genes was done with identical pipeline. (C) Based on known editing sites in representative species (D. melanogaster and C. chinensis), a combination method was proposed to identify potentially conserved RNA editing sites in target species (Apis mellifera). (D) Phylogeny of the representative species used in this study. The branch length is not proportional to the divergence. The fractions of nucleotides at orthologous positions were displayed for each species. D. melanogaster contributes most of the known RNA editing sites and thus has the highest fraction of adenosines. Gaps generally include the gaps in orthologous genes (purple) or the situation where no orthologous genes were found at all (grey).

Although the DNA-resequencing data from the matched honeybee individuals were not all available, we argue that if an A-to-G variation was observed in the honeybee transcriptome, this variation is more likely to be a conserved RNA editing site rather than SNP because it is very rare to see orthologous sites to have DNA polymorphism in one species (honeybee) and RNA editing in another species (fly). A more plausible explanation is that this conserved editing site was inherited from the common ancestor of those extant species (like in Figure 1C if a candidate adenosine in honeybee is edited).

We collected the 3308 CDS editing sites in Drosophila melanogaster (Diptera) and 113 CDS editing sites in Coridius chinensis (Hemiptera) (see Materials and Methods). According to the insect phylogeny, C. chinensis is an outgroup of bees and Drosophila (Figures 1C,D). Since we are going to compare the honeybee editing sites identified from the new method with the previous version, we should not include our de novo identified honeybee sites into the candidate list (or the old version is definitely a subset of the new editing sites). Totally a union of 3699 editing sites in CDS were obtained. Although these 3699 sites are edited in at least one species of D. melanogaster and C. chinensis, the orthologous genome positions in other species might contain non-adenosine nucleotides or even gaps (Figure 1D). For this union of 3699 sites, we calculated the fraction of each nucleotide among the orthologous positions of D. melanogaster, C. chinensis, A. mellifera, and B. terrestris. Particularly, in honeybee, we obtained 1548 adenosines (41.8%), meaning that 41.8% mapped to A in the honeybee genome. Similarly, 234 cytidines (6.3%), 355 guanosines (9.6%), 271 thymidines (7.3%), and 1291 gaps (34.9%) were obtained. The gaps include 678 sites without orthologous genes in honeybee, and 613 sites deleted in the orthologous genes in honeybee (Figure 1D). All these information and statistics can be extracted from the CDS alignments (Figures 1B,C). The 1548 adenosines in honeybee genome were defined as ‘candidate adenosine sites’. However, the existence of genomic adenosines does not ensure the occurrence of RNA editing events. We therefore need to map the transcriptomic data to the corresponding sites to detect variations in RNAs.

Mapping the honeybee transcriptome to the candidate editing positions in CDS

To determine whether the conserved genomic adenosines undergo A-to-I RNA editing, we aligned the RNA-Seq reads to the whole CDS sequences of honeybee and extracted the information on candidate adenosine sites (for the rationale of the whole procedure, please refer to Materials and Methods). The A-to-G variations in RNAs were the potentially conserved A-to-I editing sites between honeybee and at least one of the three representative species. Given the prior knowledge that these adenosines are already edited in other species, then this orthology-based methodology is able to retrieve the conserved editing sites that have been discarded by the stringent traditional de novo identification procedure.

Using REDItools to treat the transcriptome data, the sequencing coverage (Cov), alternative reads count (Alt), and editing level (L) were recorded for each of the candidate adenosine sites. To avoid false positive events by sequencing errors, we performed a binomial test on variation sites in each sample followed by multiple testing correction [23,37] (Materials and Methods). The A-to-G variations with FDR < 0.05 in any of the samples were regarded as reliable A-to-I RNA editing sites in honeybee. Then, RNA editing sites from different samples were combined and analysed.

Novel conserved recoding sites between honeybee and other insects

In total, 10.2 M to 139.2 M reads in the honeybee samples were successfully aligned to the reference CDS sequences (Table 2). Editing signal was detected in 88 adenosine sites in at least one sample of honeybee, including 72 recoding sites and 16 synonymous sites (Supplementary Table S1). The numbers of editing sites in a single sample ranged from 2 to 51 (2 to 40 for recoding sites and 0 to 12 for synonymous sites) (Figure 2A and Table 2). Due to the much larger library size of worker samples compared to drone samples, more conserved editing sites were identified in workers. Analogous to previous observations in insects, recoding sites are much more abundant than the synonymous sites among the inter-species conserved editing sites [38]. The average editing levels were 0.10 and 0.076 for nonsynonymous sites in workers and drones, and 0.079 and 0.075 for synonymous sites in workers and drones, respectively (Supplementary Table S1). Note that in D. melanogaster we previously found one editing site at stop-codon [39], changing TAG to TGG (Trp) and leading to stop-codon read-through. But this site failed to find a precise orthologous site in honeybee genome, preventing us from knowing its conservation in honeybee. Nevertheless, for all the annotated CDSs in honeybee, we investigated whether RNA editing events were found at stop-codons (Materials and Methods). A-to-I(G) mutations at TAG or TGA abolish the stop-codon while a single A-to-I(G) mutation at TAA maintains the stop-codon. However, no editing signal was observed at stop-codons in honeybee, suggesting that the few occurrences of such editing in other insects are generally non-conserved. Figure 2. Identification of conserved RNA editing sites in honeybee. (A) The numbers of conserved editing sites in each honeybee sample based on the orthologous site methodology. Recoding and synonymous sites were shown separately. (B) Conserved editing sites between honeybee and D. melanogaster. The known sites identified by previous study and the new sites identified in this study were shown separately. (C) The conservation of four recoding events in gene CG14616 (lethal-1 G0196) of two fruitflies (D. melanogaster and D. simulans) and two bees (A. mellifera and B. terrestris). The phylogeny is unscaled. The codons and AA changes were displayed for each site. ‘…’ means the codons are away from each other and ‘–’ means the codons are adjacent to each other. The averaged editing levels were obtained from brains of Drosophila adults [39] and worker bees. The IGV screenshots of representative honeybee and bumblebee samples were provided next to the diagram. The locations of the regions were indicated by CDS ID and coordinates. The editing levels were reported by REDItools with reads with mapping quality (q) ≥20 and bases with quality (Q) ≥30. Sanger sequencing was performed to validate the four RNA editing sites in honeybee gene XM_026444942 (fly ortholog CG14616). DNA and RNA (cDNA) were sequenced for each site. The Sanger traces and RNA editing levels were shown.

Table 2. Library status and conserved editing site of each honeybee sample.

Sample	Reads mapped to CDS	Nonsyn editing sites	Syn editing sites	
Drone head	10,238,984	9	2	
Drone thorax	12,156,981	4	1	
Drone abdomen	11,398,360	2	0	
Nurse brain 1	81,341,105	29	6	
Nurse brain 2	139,237,163	39	12	
Nurse brain 3	64,328,950	29	9	
Nurse brain 4	61,848,130	32	5	
Nurse brain 5	60,436,208	31	8	
Nurse brain 6	104,079,958	36	8	
Forager brain 1	67,403,058	36	5	
Forager brain 2	85,587,056	40	11	
Forager brain 3	45,810,910	33	8	
Forager brain 4	67,231,732	37	6	
Forager brain 5	77,427,424	35	10	
Forager brain 6	56,491,314	33	6	

Among the 88 edited adenosines detected in honeybee, 82 sites (67 recoding and 15 synonymous) were conserved between honeybee and D. melanogaster (meaning that the 88–82 = 6 remaining sites were conserved between honeybee and C. chinensis). Compared to the previous version of honeybee editing sites where only five sites were conserved between honeybee and D. melanogaster (two recoding and one synonymous site in gene Shab and two recoding sites in gene qvr) [34], our current work has identified 77 novel conserved editing sites, including 63 recoding sites and 14 synonymous sites (Supplementary Figure S1 and Figure 2B). These novel cases highlight the success and efficiency of our new methodology in finding previously unnoted conserved editing sites.

Intriguingly, one may raise a question that the de novo identified editing sites can also be focused on the conserved CDS sites, then what is the point of this orthology-based pipeline? In fact, in the traditional pipeline for de novo identification of RNA editing sites (Supplementary Figure S1), all the 82 honeybee-D. melanogaster conserved editing sites have been detected in the initial step (calling RNA-Seq variants against the reference genome). However, due to the low A>G% fraction at the initial step, those A-to-G variants could not be automatically regarded as RNA editing sites [24,26–30]. Then, it is obligated to perform additional filters to enrich the A-to-G variants. The stringent criteria lead to the dramatic decrease of total variants and the loss of a number of true-positive A-to-I RNA editing sites (Suplementary Figure S1). But our orthology-based methodology nicely complements the de novo identification. This method is based on prior knowledge that those particular sites are already edited in D. melanogaster, then the A-to-G variants at the orthologous sites in honeybee transcriptomes are unlikely to be artefacts, and finally these RNA editing sites are retrieved.

Moreover, even based on the 82 honeybee-fruitfly conserved editing sites, 22 sites (19 recoding and three synonymous) were also edited in bumblebee. This number already exceeded the previously reported nine conserved editing sites between honeybee and bumblebee [34], let alone only five sites were previously found to be conserved across fruitfly and two bees. Intriguingly, the analysis in bumblebee seems irrelevant to our identification of RNA editing sites in honeybee. However, regarding the confounding factors that most RNA editing studies face (the false-positive variants caused by sequencing errors, SNPs, or mis-alignments) [23], a smart way to erase the doubt is to show that the RNA variants appear in multiple species. Concurrence of sequencing errors or SNPs across multiple species/samples should be very rare, and therefore the conservation of ‘RNA variations’ in bumblebee increases the reliability of the novel RNA editing sites found in honeybee.

For example, gene CG14616 (lethal-1 G0196) encode a bifunctional inositol kinase that regulates apoptosis, vesicle trafficking, cytoskeletal dynamics, and exocytosis [40]. This gene has as many as four highly conserved recoding sites across Drosophila, honeybee, and bumblebee (Figure 2C). The 2nd ~ 4th recoding sites had identical types of AA changes in all species, while the 1st recoding site had different pre-editing AAs between fruitflies and bees but the post-edited AA was the same, a terminology called ‘conserved editing with non-conserved recoding’ [41,42]. The combinations of the four recoding sites largely increased the isoform diversity of this kinase. The 2nd Ser>Gly recoding site had editing level >0.5 in all four species, suggesting a potential essentiality for keeping this site editable during evolution. The post-edited kinase isoform might have particular functions under distinct conditions. We re-emphasize that although the matched DNA resequencing data were not available for all species, this striking conservation pattern was unlikely caused by artefacts like SNPs or sequencing errors because (1) balancing selection on the same heterozygous SNP in all four species should be extremely rare, and the heterozygous SNPs should have level ≈ 0.5 in the RNA-Seq of all species rather than in only one or two species; (2) sequencing errors should appear on all nucleotides rather than exclusively introducing A-to-G substitutions in the RNA-Seq (Figure 2C).

Moreover, Sanger sequencing was performed to validate the four RNA editing sites in honeybee gene XM_026444942 (fly ortholog CG14616). New honeybees were collected and both DNA and RNA (cDNA) from bee heads were sequenced (Materials and Methods). We could see that for the most highly edited site2 (Ser>Gly), RNA has both A and G alleles while DNA only has A; and the RNA editing level in Sanger sequencing is 0.56, very close to the level obtained in NGS (Figure 2C). This suggests the reliability of our pipeline in finding conserved RNA editing sites in CDS. For the other three lowly edited sites with NGS editing level <0.03, Sanger sequencing is difficult to precisely quantify the editing level due to the background noise (see Materials and Methods for a more detailed explanation), but the estimated editing levels from the Sanger traces are generally similar to the level in NGS (Figure 2C). Despite the limitation in the validation of lowly edited sites, we propose that the future functional studies should prioritize the highly edited sites conserved across multiple species (like the site2 Ser>Gly recoding) since they are more likely to be functional and essential.

Apart from finding conserved editing sites, another use of our orthology-based methodology is to identify species-specific editing sites. We acknowledge that conserved editing sites might be intuitively thought to be more functional than non-conserved ones, but our approach indeed could find both categories. Among the 3308 known RNA editing sites in CDS of D. melanogaster, 1524 sites matched to an adenosine in the honeybee genome and thus the remaining 1784 editing sites are specific to D. melanogaster. Then, among the 1524 adenosines in honeybee, 549 sites have RNA-Seq coverage ≥100 and no editing signal detected. This suggests that the 1784 + 549 = 2333 sites are likely the D. melanogaster-specific sites.

Caste-specific recoding sites in honeybee among the conserved CDS sites

One of the most amazing features of eusocial insects is to generate behaviourally distinct individuals from the same genome sequence. This phenomenon highlights the importance of potential regulatory mechanisms beyond the genome. Particularly, what drives the differentiation of different castes and sub-castes? It is possible that RNA editing might be an essential post-transcriptional approach to shape the behaviour and phenotype of honeybees. Caste-specific RNA editing is an interesting issue to be investigated. Indeed, the original honeybee RNA editing paper [34] and the previous bumblebee RNA editing paper [33] have done such differential editing analysis in the corresponding species. Here, we would like to use the particular examples of caste-specific editing to show that the RNA editing events found by orthology-based method are reliable: the caste-specific editing sites with differential editing levels between two castes are the best evidence to show that the variation detected in RNA-Seq is unlikely caused by heterozygous SNPs (because heterozygous SNPs should have similar allele frequencies ≈ 0.5 in different samples) or sequencing errors (very low frequency in all samples, if any).

Among the total 88 conserved editing sites in honeybee, we set out to identify differential editing sites (DES) between different castes (drone versus worker) or sub-castes (nurse versus forager) (see Materials and Methods for details). We identified 2 DESs between drone and workers where drone had higher editing levels, and both sites were recoding sites. One AGC>GGC (Ser>Gly) recoding site was located in gene RhoBTB (Rho-related BTB domain) which encodes a protein involved in anti-parasitoid immune response. Given the biological nature of drones and workers, it is possible that the two castes are subjected to different levels of threats from the environmental parasites, and then the recoding on RhoBTB gene might be programed to defend the host against the parasites. Another ACA>GCA (Thr>Ala) recoding DES was located in gene Cul4 (Cullin 4) encoding a molecular scaffold for the CRL4 E3 ubiquitin ligase complex.

Then, five DESs were found between nurses and foragers, and again, all of them were recoding sites. One AGT>GGT (Ser>Gly) recoding site was located in a well-known neuronal gene qvr (quiver) which encodes a Ly-6 protein that modulates the trafficking and activity of membrane protein targets. Since the difference between nurses and foragers is mainly reflected by the behaviour rather than morphology, the differential recoding on neuronal-related genes might nicely explain the phenotypic divergence. Next, we will highlight a special DES of our interest and infer the biological significance behind caste-specific RNA editing.

Caste-specific auto-recoding in honeybee Adar gene and the nearby passenger editing events

Adar is a typical example of auto feedback regulatory enzyme. The best-known mechanism was in Drosophila where a Ser>Gly auto-recoding site in deaminase domain reduced the catalytic activity [43,44]. Then, an ATA>ATG (Ile>Met) auto-recoding site in Adar deaminase domain was found in bumblebee where the auto-recoding level was correlated with the global editing efficiency [33].

Since eusocial insects are able to create great phenotypic plasticity, RNA editing mechanism might be an essential approach that should be precisely controlled. Indeed, we found that the Adar Ile>Met auto-recoding level was significantly higher in drones compared to workers, and also significantly higher in foragers compared to nurses (Figure 3A). This caste-specificity highlights the possibility that RNA editing is well suited to contribute to the transcriptomic and phenotypic plasticity of eusocial honeybees. Figure 3. Auto-editing sites in honeybee adar gene. (A) Differential editing levels observed at Ile>Met recoding site. The difference between drone head and worker brains was judged by Fisher’s exact test. The difference between six nurses and six foragers was determined by T-test. (B) The sequence context of Adar Ile>Met auto-recoding site. Six auto-editing sites and the functional consequence were Illustrated. Recoding sites are in red and a synonymous site is in blue. (C) IGV visualization of the passenger editing events linked to the Ile>Met recoding site in Adar CDS. Drone head and nurse brain 1 were shown as examples. The linkage disequilibrium (LD, r2 and p value) between site #2 and site #6 was demonstrated. To be consistent with the REDItools output, here, only the reads with mapping quality (q) ≥20 and bases with quality (Q) ≥30 were displayed and used for calculating LD.

Moreover, we identified five novel auto-editing sites in Adar, four are recoding sites and one is a synonymous site (Figure 3B). The Ile>Met site and its two upstream codons made up a 9 nt sequence GAC-AAA-ATA, RNA editing was observed at GAC (site #1, recoding, Asp>Gly), AAA (site #2, recoding, Lys>Glu), AAA (site #3, recoding, Lys>Arg), AAA (site #4, synonymous, Lys), and ATA (site #5, recoding, Ile>Val; or Met>Val if the 3rd codon position was already edited). The Ile>Met site itself would be site #6 sequentially (Figure 3B).

However, only the Ile>Met recoding site at ATA constantly appeared in all the tested samples while the other auto-editing sites only randomly occurred in a few samples (Figure 3C and Supplementary Figure S2). This raises a possibility that Ile>Met recoding is the main target of Adar and the other nearby editing sites are passengers or byproducts [45]. This assumption is further supported by the observation of strong linkage between the Ile>Met site and nearby editing sites (Figure 3C). In sample nurse brain 1, three adenosines upstream the Ile>Met site were edited (not edited in the same reads), two of which were nonsynonymous and one was synonymous editing. Interestingly, all the upstream editing events were linked to the editing events at Ile>Met site; in other words, the upstream editing sites is only edited when the Ile>Met site is edited (Figure 3C). This ‘single direction complete linkage’ was also widely observed in other samples without any exceptions (Supplementary Figure S2), suggesting a passenger role of the additional editing sites. For example, in drone head, a significant and strong linkage disequilibrium (LD, r2 = 0.64, p = 6.9E–4) [46] was observed between site #2 and the Ile>Met site (Figure 3C).

These results reveal a tolerance of promiscuous editing by Adar. The Asp>Gly and Lys>Glu recoding events upstream the Ile>Met site might cause undesired AA substitutions to the Adar protein, and the ATA editing could directly change the AA to Val no matter whether the 3rd codon position is edited (Figure 3B). We surmise that the benefit of Ile>Met editing exceeds the disadvantage of promiscuous editing at nearby regions, so that natural selection has maintained this mechanism where Ile>Met recoding is constantly highly edited in all samples while the nearby passenger editing events were lowly edited varied randomly across different samples (Figure 3C and Supplementary Figure S2). Nevertheless, we do not rule out an alternative explanation that this seemingly random passenger editing events are actually contributing to the great plasticity and inter-individual differences of eusocial insects.

Caste-specific Ile>Met recoding and the passenger editing events are conserved in bumblebee

In addition to the interesting auto-recoding in honeybee, we obtained a more striking result when we looked at the Ile>Met auto-editing in bumblebee Adar gene (Figure 4). In the 24 brain samples of bumblebee (eight nurses, eight foragers, and eight newly emerged callow workers) [33], we detected exactly the six editing sites as we found in honeybee Adar. This observation directly adds five new conserved editing sites between honeybee and bumblebee (except the Ile>Met recoding is already a known conserved site). The genome sequence of the three consecutive codons was identical between honeybee and bumblebee where site #6 is the canonical Ile>Met (ATA>ATG) recoding sites with the highest editing level (Figure 4 and Supplementary Figure S3). Interestingly, in bumblebee, site #2 (AAA>GAA, Lys>Glu) also had significant linkage (p = 0.0015) with the Ile>Met site (Figure 4A), suggesting that not only the editing sites but also the potential significance of this linkage phenomenon were highly conserved between the two bees. Moreover, through manual inspection, we confirmed that the occurrences of site #1 to site #5 were strictly dependent on (linked to) the Ile>Met editing, meaning that there were no RNA reads with site #1 to site #5 edited but with site #6 unedited (Figure 4A). Since this intriguing pattern was also conserved between two bees, it is encouraging to believe that the Ile>Met recoding is the driver site and the nearby editing events are passengers. Figure 4. Conservation of adar auto-editing and caste-specificity between bumblebee and honeybee. (A) IGV visualization of the six auto-editing sites in Adar CDS of bumblebee. Site #6 is the Ile>Met (ATA>ATG) recoding site. The linkage disequilibrium (LD, r2 and p value) between site #2 and site #6 was demonstrated for sample nurse 1. Not all reads were shown due to space limitation. Then, only the reads with mapping quality (q) ≥20 and bases with quality (Q) ≥30 were used for calculating LD. The IGV screenshots for other worker samples were shown in Supplementary Figure S3. (B) Differential editing level of the Ile>Met recoding site between different bumblebee workers. Callow means newly emerged callow workers. p values were obtained by one-tailed T-tests.

Then, we questioned whether the Ile>Met recoding is caste-specific in bumblebee (Figure 4B). Using the 24 samples from three different sub-castes of workers, we found that the newly emerged callow workers and the differentiated foragers had similar editing levels (median between 0.40 ~ 0.42) while the nurses had significantly higher editing levels (Figure 4B). Since the callow workers were used as a control, it could be inferred that the foragers have maintained an original state at the Ile>Met recoding site, while nurses have elevated the Ile>Met editing level during caste differentiation. Note that although the caste-specificity of the Ile>Met recoding level is conserved between honeybee and bumblebee, their directions are opposite: in honeybees, foragers have higher Ile>Met recoding levels. This highlights the flexibility of RNA editing and its advantage in shaping the transcriptomic plasticity of eusocial insects.

The orthology-based methodology is not affected by differentially expressed genes

Similar to the traditional de novo identification of RNA editing sites, detection bias is a potential concern in virtually all RNA editing or gene expressional analyses, including our orthology-based editing detection. To understand how the gene expression pattern across different sample might affect the identification of conserved editing sites, we performed differential expression analysis between six nurses and six foragers. We found 873 gene with FDR < 0.05 and |log2foldchange| > 0.5 among the total 23,467 genes (Materials and Methods). We therefore obtained 873 differentially expressed genes (DEG) and 22,564 non-DEG. The totally 88 editing sites detected by our orthology-based methodology belong to 69 unique genes, among which 2 genes are DEG (2/873 = 0.229%) and 67 genes are non-DEG (67/22564 = 0.297%). The two fractions are not significantly different (p = 0.96 under Chi-square test). This suggests that no matter how variable a gene is across different samples, the conserved editing sites are likely to be detected, at least this is true for the 88 conserved editing sites identified by our orthology-based pipeline. This observation adds confidence to the reliability of the novel RNA editing sites we found. Moreover, our Sanger validation on the particular RNA editing sites further supports their authenticity (Figure 2C). In future studies on the function and adaptation of RNA editing sites, we suggest that the conserved recoding sites with relatively high editing level should be prioritized.

Discussion

Summary of main findings

In this work, we proposed an orthologous site methodology to retrieve the previously missed conserved editing sites. Those true-positive editing sites might be discarded during the highly stringent de novo identification of RNA editing sites (which requires a high A>G% fraction among all variations), but our method benefits from a prior knowledge that a position is edited in other species and therefore we could pick up a number of editing site with less stringent filtering steps. We found a considerable amount of novel recoding sites in honeybees that were conserved across multiple insect species, providing potential candidates for future functional studies. Then, caste-specific editing sites were identified, including the interesting Ile>Met auto-recoding site in Adar CDS. The caste-specificity of this Ile>Met recoding sites, together with several nearby passenger editing sites linked to it, were conserved between honeybee and bumblebee. This conservation pattern suggests a putative role of this Adar auto-regulatory editing site in shaping the phenotypic plasticity of eusocial insects.

Further application of the orthology-based methodology

We acknowledge that in this work, no software or algorithms were new, but we are the first to use an orthology-based pipeline to identify many previously unnoticed RNA editing sites in a species. Given that the species with RNA editing studies only represent a corner of the iceberg, none of previous studies have used this method to retrieve the previously ignored RNA editing sites. We believe that this is the novelty of our work. Particularly, since honeybee released a better genome version recently (suggesting that previous RNA editing study based on old genome version might missed many true positive sites due to mapping issues), then our methodology can largely benefit the identification of new RNA editing sites in honeybees. On the other hand, we are fully aware that the aim of our current study is not the de novo identification of RNA editome in honeybee, but to use a new methodology to retrieve the conserved editing sites. De novo identification of RNA editome is continuously being done in new species [1], and this allows us to combine more and more sites in different species to replenish the union of candidate RNA editing sites. Nevertheless, we should stress that many of the de novo identification works were done in mammals [47], molluscs [15], and early-diverging metazoans [48], while RNA editing in the insect clade was less investigated [41] although the number of insect species makes up more than half of the animals. Intriguingly, insect RNA recoding sites were thought to have a relatively higher conservation level compared to the recoding sites in mammals. For example, roughly two thirds of the recoding sites in D. melanogaster were detected in other Drosophila species [38,39,49], but the mammalian recoding sites were poorly conserved given the large basal number of RNA sites in each mammal [21,50]. Therefore, we anticipate that when a balanced attention is paid to the RNA editing in insect species, the list of candidate editing sites could be dramatically expanded. This not only reveals some previously unheeded editing events but also allow the discovery of more highly conserved and functional recoding sites for future experimental validation.

Manual inspection reveals novel editing sites independent of any detection strategies

The orthology-based methodology aims at finding conserved editing sites in CDS, but might ignore species-specific and sample-specific editing sites. Then, the manual inspection via IGV or sequence alignment file is helpful in finding such lowly conserved and lowly edited editing events. The several additional auto-editing sites found in honeybee Adar gene are nice examples of successful manual inspection. Notably, one may argue that the hyper-editing pipeline should be able to discover these weak editing events [51]. However, (1) Typical hyper-editing events occurs in lowly expressed regions and the region almost only expresses the edited molecules [9]. But Adar in brains is highly expressed and the auto-edited molecules only make up a small fraction of the total reads as reflected by the low editing level. (2) The identification of hyper-editing usually requires multiple clustered editing events within the same molecule. But for the Adar auto-editing sites in our cases, although there are multiple positions being editable, one RNA read usually only contains at most two editing events (one Ile>Met editing event linked to another nearby editing event). Computationally, this situation might not reach the hyper-editing threshold. As a consequence, some ‘lowly conserved and lowly edited sites in highly expressed genes’ can only be discovered by manual inspection.

What can we learn from the passenger editing events?

Having said that, those editing sites found by manual inspection might be less functional than the conserved recoding sites found by the orthology-based methodology. But we are quite certain that those passenger editing events around Ile>Met site were not sequencing errors because the errors should be randomly distributed among the reads and will not show strong linkage to the Ile>Met editing events. It is worth thinking why and how the organisms can tolerate these slightly deleterious byproducts. For example, the Asp>Gly and Lys>Glu editing events linked to the Adar Ile>Met recoding makes another AA change to the AdarMet protein, and the Met>Val editing at the same codon of the Ile>Met recoding site directly changes the AdarMet isoform to a functionally unknown AdarVal isoform. These undesired alternations are unfavourable, and question comes that why should organisms take this risk? We propose the following explanation. It should be noted that although the passenger editing events are always linked to the Ile>Met editing event, the Ile>Met editing event is not always accompanied by a passenger editing event. That is to say, most of the Ile>Met edited Adar mRNAs (~25% of the total expressed Adar mRNA) faithfully produce an AdarMet isoform and only a few produce an abnormal one. The benefit of the high level Ile>Met recoding (e.g. controlling the global editing efficiency) far exceeds the damage caused by the few passenger editing events (e.g. the waste of some energy/resource and the production of some potentially toxic proteins [52]). Therefore, natural selection does not have to suppress the occurrence of such byproducts.

Summary and emphasis

Taken together, we proposed a complementary approach to the traditional pipeline and retrieved several previously unnoticed CDS editing sites. We stress that we do not claim to develop a better methodology that comprehensively identifies all the A-to-I RNA editing sites in a given species. Neither do we claim that our orthology-based methodology can identify more editing sites than the traditional de novo identification pipeline. Instead, our key point is to show the complementarity and cooperation between the two approaches. In conclusion, from both technical and biological aspects, our works facilitate future researches on finding the functional editing sites and advance our understanding on the connection between RNA editing and the great phenotypic diversity of organisms.

Materials and methods

Data availability

The honeybee genome HAv3.1 was submitted by Wallberg et al. [53] and we downloaded the reference sequence, gene annotation, and repeat annotation from NCBI (https://www.ncbi.nlm.nih.gov/). The genome sequence and annotation of Drosophila melanogaster were downloaded from FlyBase (https://flybase.org/) version dm6.04. The genome assembly and list of A-to-I RNA editing sites of Coridius chinensis were retrieved from our previous study [41]. The known editing sites of D. melanogaster were combined from a large set of different studies [22,34,39,45,49,54–58]. A union of 7422 candidate RNA editing sites in D. melanogaster was obtained, including 3308 (44.6%) sites located in CDS, followed by 2247 (30.3%) intronic sites, 946 (12.7%) sites in UTRs, and other 921 (12.4%) sites in non-coding RNAs or unannotated genomic regions. We downloaded the transcriptomes of multiple tissues and castes of Apis mellifera with the following accession IDs: the transcriptomes of a single male drone were downloaded from Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa/). RNA-Seq for individual#1: CRX082741 for head, CRX082742 for thorax, and CRX082740 for abdomen. RNA-Seq for individual#2: CRX082744 for head, CRX082745 for thorax, and CRX082743 for abdomen. DNA-resequencing data for individuals #1 and #2 are CRR106402 and CRR106403, respectively. The brain transcriptomes of female workers were downloaded from NCBI (SRR445999-SRR446004 for six nurses and SRR446005-SRR446010 for six foragers). The drone data were single-ended 50 bp and the worker data were pair-ended 100 bp. The layout and length of the reads mainly affected the hyper-editing detection but had little effect on identifying conserved editing sites based on the orthologous site methodology. The transcriptomes of bumblebee Bombus terrestris were retrieved from a previous study [33]. Sanger sequencing raw files were provided in Supplementary Data 1.

The rationale for our orthology-based method to identify conserved CDS editing sites

Mapping RNA-Seq reads to the honeybee reference genome belongs to the traditional pipeline, aiming to identify the overall RNA editing sites in a species regardless of whether this editing is conserved across species. The traditional method has been used in the previous honeybee [34], bumblebee [33], and leaf-cutting ant [25] studies. The limitation of the traditional pipeline has been mentioned in Introduction. In contrast, the basic idea of our proposed methodology is to directly look at the sites in a species (e.g. honeybee) which corresponds to known editing sites in another species (e.g. D. melanogaster) (Figures 1B,C). The two methods are complementary and we do not claim which is a better one because they focus on different aspects, but they do share part of editing sites as we have revealed. To perform the orthology-based method, the alignment is the first step. Note that the matching of genomic locations between different species is more accurate in CDS compared to non-coding regions (because CDS alignments are more robust and more applicable/feasible compared to whole-genome alignments). Thus, we only focus on CDS editing sites.

We first aligned the CDS sequences of different species to determine the CDS coordinates in honeybee that corresponds to the known CDS editing sites in D. melanogaster (Figure 1C). Then, we only need to map the RNA-Seq reads to the honeybee CDS and see if the candidate sites are edited. Mapping RNA-Seq reads to CDS avoids the complicated issue of dealing with reads spanning splicing junctions [59] and should be more accurate than mapping RNA-Seq reads to the whole genome. But due to the trade-off between quantity and quality, most common studies still map RNA-Seq to the whole genome to cover the non-coding regions such as introns or non-coding RNAs. We only focus on potential conserved CDS editing sites, and this is why the RNA-Seq mapping step took place after the CDS alignment step. To achieve this purpose, we do not need to map the RNA-Seq reads to the whole genome. But this is only for the orthology-based methods. For the hyper-editing method, we still mapped the RNA-Seq reads to the whole genome of honeybee.

Sequence alignment, phylogeny, and candidate adenosine sites in honeybee

The CDS sequences were translated into protein and aligned using blastp v2.2.28 [60]. Default parameters were used. A complete schematic flow chart for this part is provided in Figure 1B. For example, when searching D. melanogaster-honeybee orthologous genes, we first obtained that the 3308 CDS editing sites in D. melanogaster belong to 1547 unique D. melanogaster genes (proteins). We aligned these 1547 genes (protein sequences) to the honeybee reference protein sequences, finding that 972 genes are reciprocal best. ‘Reciprocal best’ means that a honeybee gene B1 is the best hit of fly gene A1, and A1 is also the best hit of B1. Indeed, one honeybee gene might find many hits in fly, and vice versa, but we only kept the ‘reciprocal best pairs’ in this round. Blastp has a default parameter of E-value <1E–6 (lower E-value means better match) and the multiple hits are ranked by increasing E-value. We acknowledge that the multiple hits are especially prevalent within paralogs or gene families, but one can always rank the candidates with the default order in blastp and then the best match can be chosen. Conceivably, the reciprocal best pairs are high-confidence orthologous genes between honeybee and fly. However, we also consider that the reciprocal best strategy might be over stringent so we carried out the second round of blastp for the 1547–972 = 575 D. melanogaster genes remaining. We aligned these 575 genes (protein sequences) to the honeybee reference protein sequences (with the previous 972 honeybee genes excluded). We required E-value <1E–6 and matching length > 50% on both sides and picked up the best matched honeybee gene if any. By this second round blastp, we further retrieved 136 orthologous genes between fruitfly and honeybee. This still face the situation that one gene from fruitfly may match two or more genes in honeybee, but as we have clarified above, we selected the best match. The second round does not require reciprocal best because all reciprocal best pairs have already been selected in the first round. Therefore, among the 1547 edited genes in fruitfly, 972 + 136 = 1108 (71.6%) genes have found an orthologous gene in honeybee. In other words, 1108 pairs of fly-bee genes were obtained. Then the CDS and protein sequences between fruitfly and honeybee were converted into alignment format using MAFFT v7.158b [61]. Notably, to increase accuracy, CDSs were aligned referring to the pre-aligned protein sequences. Finally, the coordinates of known CDS editing sites in D. melanogaster were projected to the honeybee CDS based on the alignment file. The mapping of C. chinensis known editing sites to the honeybee CDS was done by the same strategy.

To show how well the alignments are, we take the 1108 D. melanogaster-honeybee orthologous genes as an example. For each pair of gene, we first calculated the ‘proportion of alignment without gap’. This is a direct measurement of how well the alignment is. Higher proportion indicates higher quality and reliability of the alignment. We found that the 1108 pairs of orthologous genes had a median value of 78.8% (Supplementary Figure S4, mean = 75.7%, S.E. = 0.5%). Then, for the ungapped region, the median identity is 52.0% (Supplementary Figure S4, mean = 53.3%, S.E. = 0.4%). The alignments are already of high quality given the great divergence time between Diptera and Hymenoptera. The high proportions and identities are likely due to the fact that the edited genes in Drosophila are generally of high conservation level [38,49] so that the orthologous genes in other insects are likely to have high similarity.

According to the sequence alignment, we made a union set of 3699 candidate RNA editing sites from D. melanogaster and C. chinensis. These 3699 orthologous positions were projected to the CDS coordinate of honeybee. Among them, 1548 (41.8%) sites were adenosines in honeybee, representing the ‘candidate adenosine sites’ of our study. Notably, there are other two Hymenoptera species leaf-cutting ant Acromyrmex echinatior [25] and bumblebee Bombus terrestris [33] that have RNA editing sites reported. The reasons for not incorporating these candidate sites are that (1) Both species have an updated genome version. Since part of our goal is to show the power of using an improved genome assembly, we did not include these two species; (2) The bumblebee RNA-Seq data were used as a control to verify if the honeybee-fruitfly conserved editing sites were also detected in bumblebee; (3) The ant paper only provided the list of recoding sites but not the synonymous sites [25], and for the bumblebee paper, most of the provided editing sites were in non-coding regions [33], which does not add too much to our analysis.

Note that the reference CDS sequences from different species were first used for the alignment. And the alignment file (containing gaps) was used for aligning the orthologous sites and then to infer the relative positions of those sites on each ungapped CDS. In the next transcriptome mapping step, the RNA-Seq reads were still mapped to the reference CDS sequences (without gaps).

Mapping the transcriptome to the candidate editing site in CDS

We selected the brain RNA-Seq of six nurses and six foragers, and the RNA-Seq of head/thorax/abdomen of drone individual#1 were also used. The drone individuals #1 and #2 were subjected to polyA+ library construction. But the quality of individual#2 was relatively low (with <60% uniquely mapping rate to the genome), and thus we only used individual #2 for the de novo identification as mentioned below (this is to be consistent with the scheme of our previous honeybee study). In this orthology-based methodology, we only used individual#1.

We first mapped the transcriptome reads to the honeybee reference genome HAv3.1 using STAR version 2.7.6a [59]. Default parameters were used. Based on the ‘NH:i:1’ tag in the alignment file, we extracted the reads that have mapped to a single location in the genome, termed ‘unique mappers’. The number of unique mappers of each library was regarded as library size.

We then mapped the unique mappers to the reference CDS sequence using BWA mem version 0.7.17-r1198 [62]. Default parameters were used. The numbers of reads successfully aligned to CDS would be informative and were recorded in Table 2. The sequencing coverage (Cov) and alternative reads count (Alt) on candidate adenosine sites were extracted using REDItools v2.0 [63] with default parameter which requires mapping quality (q) ≥ 10 and base quality (Q) ≥ 30. If a variation was found at an candidate adenosine site by REDItools, we would confirm the reliability of the variation by GATK HaplotypeCaller version 4.3.0.0 [64]. Only the variations supported by both REDItools and GATK were maintained. For A-to-G variations on candidate adenosine sites, editing level was defined as G/(G+A). We further filtered the sites by a binomial test [23] on the Cov and Alt values and required FDR < 0.05 after multiple testing correction [37]. The A-to-I RNA editing sites passing the filter in any of the samples were regarded as an authentic editing site in honeybee. The bumblebee transcriptome and RNA editing sites were treated with identical approaches. The orthology of the sites were known from the CDS alignment in the previous steps.

Traditional pipeline for identification of A-to-I RNA editing sites

The traditional pipeline for identification of A-to-I RNA editing sites, namely the de novo identification, has been well described by different studies including our owns [23,35,65]. Please also refer to our previous honeybee paper to understand how the data were produced. Here, we used the six RNA-Seq from two drone individuals that were subjected to polyA+ library construction. For data analysis, in brief, based on the RNA-Seq alignment mapped to the reference genome HAv3.1 using STAR version 2.7.6a [59], duplications in the alignment were first removed by Picard tool v1.124 (https://broadinstitute.github.io/picard/index.html) and then the genome-wide variations were extracted with REDItools v2.0 [63]. Then, DNA-Seq reads were mapped to the reference genome using BWA mem version 0.7.17-r1198 [62]. Similar to the orthology-based method, only the variations supported by both REDItools and GATK [64] were maintained, and the sequencing coverage (Cov), alternative reads count (Alt), editing level, binomial test, and FDR were also identically defined. In general, intergenic variations are not considered because (1) these regions should not have RNA-Seq signals in theory; (2) even we believe that some real genic regions are mis-identified as intergenic regions in the reference genome, we cannot determine the direction of the variations in RNA as we do not know which strand does the gene belong to: a T-to-C variation in the antisense strand RNA might be mis-regarded as an A-to-G variation in the sense strand RNA. But genic regions do not have this problem as we know the strand information of each gene.

A reliable RNA editing site (RDD) in a sample can be defined as having RNA-Seq variation with FDR < 0.05 and meanwhile be supported by ≥10 reads in DNA-Seq and no alternative alleles appear in DNA. However, in some cases, a single sample/individual cannot produce a satisfactorily high A>G% fraction (for unknown reasons), then one may require a variation site to appear in multiple samples/individuals that all meet this set of criteria. Then the A>G% will be largely improved. The reason for not using the public worker (nurses and foragers) data in the de novo identification is the lack of a matched DNA-resequencing library from the same individuals, preventing us from the accurate exclusion of SNPs.

Gene expression analysis

Gene expression was measured by RPKM (reads per kilobase per million mapped reads). The number of reads mapped to each gene was counted with featureCounts [66]. Only exonic reads were used. For each gene in each sample, RPKM = number of reads/length(Kb)/library size. The differential expression analysis between six nurses and six foragers was accomplished by DESeq2 v1.34.0 [67]. Default parameters were used. The genes with FDR < 0.05 and |log2foldchange| > 0.5 were regarded as differentially expressed genes (DEGs). A normalized reads count for each gene per sample was generated by the software, and the mean value for six nurses and mean value for six foragers were used to plot the expressional comparison between the two sub-castes.

Differential editing site (DES)

DESs were identified between different castes or sub-castes. For the comparison between drone head and worker brains, since only one drone head sample was available, T-test was not applicable, we therefore tested the difference between editing level in drone head and editing level in 12 pooled worker samples. The coverage (Cov) and alternative allele count (Alt, G allele) of two castes was compared by Fisher’s exact test. The p values were adjusted for multiple testing correction to obtain the false discovery rate (FDR) [37]. DESs between drone and workers were defined as editing sites with FDR < 0.05 in the comparison.

For the comparison between six nurses and six foragers (brains), a six versus six T-test was applied to the editing levels of each site. p values were adjusted for multiple testing correction [37]. Meanwhile, a p value followed by FDR was obtained from the Fisher’s exact test between pooled reads of nurses and pooled reads of foragers. DESs between nurses and foragers were defined as editing sites with FDR < 0.05 in the T-test or FDR < 0.05 in the Fisher’s exact test.

Linkage between RNA editing sites

The linkage between RNA editing sites was calculated according to the original literature that introduced the formula of linkage disequilibrium (LD) [46]. For two editing sites, all we need to know is the numbers of reads supporting the four haplotypes (AA, AG, GA, and GG) to calculate the LD. The reads were manually sorted and counted in the IGV console. Then, we refined the haplotype frequencies to match the Alt and Cov counts reported by REDItools, which means, only show the reads with mapping quality (q) ≥ 20 and variation sites with base quality (Q) ≥ 30.

Sanger sequencing validation

Newly emerged female workers of honeybees (not yet differentiated into foragers or nurses) were collected from beekeepers in Beijing. To validate whether the candidate sites are edited and confirm the accuracy of the editing level, we performed Sanger sequencing on PCR-amplified genomic DNA (gDNA) and cDNA sequences. For cDNA synthesis, 500 ng of total RNA was revered transcribed using PrimeScriptTM RT reagent Kit with gDNA Eraser Kit (TaKaRa), following the manufacturer’s instructions. Primer sequences are listed in Table 3.Table 3. PCR primers used for Sanger sequencing.

Gene ID	DNA/cDNA	F/R	Primer sequences	
XM_026444942	DNA	F	CATTTTCCCGACTCCAGCCT	
XM_026444942	DNA	R	AGCTAAGACTCACCCGCAAC	
XM_026444942	cDNA	F	CGACACACGACCTTCGATCA	
XM_026444942	cDNA	R	ACCTGACAGCTTCCAAGTCG	
F: Forward primer; R: Reverse primer.

A 25 ul PCR reaction comprised EmeraldAmp® Max PCR Master Mix (TaKaRa), 100 ng of gDNA (or 5 ng of cDNA) template, and 10 μM each of forward and reverse primers. The PCR program was set as follows: 95°C for 1 min, followed by 40 cycles of 95°C for 20 s, 54°C for 30 s, and 68°C for 30 s, with a final extension at 72°C for 5 min. Primers were synthesized in Sangon Biotech (Shanghai) Co., Ltd., and Sanger sequencing was conducted by Beijing Tsingke Biotech Co., Ltd. Following our previous study [41], evaluation of RNA editing level involved measuring peak area from Sanger sequencing traces using SnapGene software (https://www.snapgene.com/). Editing level at a particular site = area of G/sum of area of A and G. However, editing levels lower than 0.05 is difficult to quantify using Sanger sequencing due to the intrinsic background noise. This is similar to the inevitable sequencing errors in NGS. The Sanger traces of nearby bases might affect the estimation on focal site. For example, for site3 in Figure 2C, the editing level in NGS is only 0.005 but the level in Sanger traces appears to be 0.016, possibly affected by the traces of the downstream guanosine. In other words, for the extremely lowly edited sites in NGS, Sanger sequencing cannot completely prove or disprove the existence of RNA editing. But for site2 with editing level around 0.5, NGS and Sanger sequencing are highly consistent.

Statistics and graphical works

Statistics and graphical works were accomplished in R language (version 3.6.3). The differential analysis was performed by DESeq2 v1.34.0 as described above. Other statistical tests or graphs were accomplished by intrinsic command line in R like ‘fisher.test’, ‘mean’, ‘plot’, ‘barplot’, ‘boxplot’.

Supplementary Material

Supplementary_FigureS.pdf

TableS1.xlsx

Acknowledgments

We thank the National Natural Science Foundation of China, the Young Elite Scientist Sponsorship Program by CAST, the Young Elite Scientist Sponsorship Program by BAST, and the 2115 Talent Development Program of China Agricultural University for the financial support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The honeybee genome HAv3.1 was submitted by Wallberg et al. [53] and we downloaded the reference sequence, gene annotation, and repeat annotation from NCBI (https://www.ncbi.nlm.nih.gov/). The genome sequence and annotation of Drosophila melanogaster were downloaded from FlyBase (https://flybase.org/) version dm6.04. The genome assembly and list of A-to-I RNA editing sites of Coridius chinensis were retrieved from our previous study [41]. The known editing sites of D. melanogaster were combined from a large set of different studies [22,34,39,45,49,54–58]. A union of 7422 candidate RNA editing sites in D. melanogaster was obtained, including 3308 (44.6%) sites located in CDS, followed by 2247 (30.3%) intronic sites, 946 (12.7%) sites in UTRs, and other 921 (12.4%) sites in non-coding RNAs or unannotated genomic regions. We downloaded the transcriptomes of multiple tissues and castes of Apis mellifera with the following accession IDs: the transcriptomes of a single male drone were downloaded from Genome Sequence Archive (GSA, https://ngdc.cncb.ac.cn/gsa/). RNA-Seq for individual#1: CRX082741 for head, CRX082742 for thorax, and CRX082740 for abdomen. RNA-Seq for individual#2: CRX082744 for head, CRX082745 for thorax, and CRX082743 for abdomen. DNA-resequencing data for individuals #1 and #2 are CRR106402 and CRR106403, respectively. The brain transcriptomes of female workers were downloaded from NCBI (SRR445999-SRR446004 for six nurses and SRR446005-SRR446010 for six foragers). The drone data were single-ended 50 bp and the worker data were pair-ended 100 bp. The layout and length of the reads mainly affected the hyper-editing detection but had little effect on identifying conserved editing sites based on the orthologous site methodology. The transcriptomes of bumblebee Bombus terrestris were retrieved from a previous study [33]. Sanger sequencing raw files were provided in Supplementary Data 1.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15476286.2024.2397757

Abbreviations

AA amino acid.

A-to-I adenosine-to-inosine.

ADAR adenosine deaminase acting on RNA.

CDS coding sequence.

DEG differentially expressed gene.

DES differential editing sites.

dsRNA double-stranded RNA.

edIle editable isoleucine codon.

edSer editable serine codon.

FDR false discovery rate.

MM mismatch.

NGS next generation sequencing.

RDD RNA-DNA difference.

RPKM reads per kilobase per million mapped reads.

S.E. standard error.

SNP single nucleotide polymorphism.

unIle uneditable isoleucine codon.

unSer uneditable serine codon.

Alt alternative allele count.

Cov coverage.

L editing level.

Authors’ contributions

Conceptualization & supervision: Y.D., W.C., and H.L.

Data analysis: Y.D., T.Z., J.L., C.Z., and L.M.

Writing – original draft: Y.D., W.C., and H.L.

Writing – review & editing: T.Z., J.L., C.Z., L.M., F.S., T.L, W.C., H.L., and Y.D.

All authors approved the submission of this manuscript.
==== Refs
References

[1] Zhang P, Zhu Y, Guo Q, et al. On the origin and evolution of RNA editing in metazoans. Cell Rep. 2023;42 (2 ):112112. doi: 10.1016/j.celrep.2023.112112 36795564
[2] Duan Y, Ma L, Song F, et al. Autorecoding A-to-I RNA editing sites in the Adar gene underwent compensatory gains and losses in major insect clades. RNA. 2023;29 (10 ):1509–1519. doi: 10.1261/rna.079682.123 37451866
[3] Ma L, Duan Y, Wu Y, et al. Comparative genomic analyses on assassin bug Rhynocoris fuscipes (Hemiptera: reduviidae) reveal genetic bases governing the diet-shift. iScience. 2024;27 (8 ):110411. doi: 10.1016/j.isci.2024.110411 39108731
[4] Bian Z, Ni Y, Xu JR, et al. A-to-I mRNA editing in fungi: occurrence, function, and evolution. Cell Mol Life Sci. 2019;76 (2 ):329–340. doi: 10.1007/s00018-018-2936-3 30302531
[5] Feng C, Xin K, Du Y, et al. Unveiling the A-to-I mRNA editing machinery and its regulation and evolution in fungi. Nat Commun. 2024;15 (1 ):3934. doi: 10.1038/s41467-024-48336-8 38729938
[6] Liao W, Nie W, Ahmad I, et al. The occurrence, characteristics, and adaptation of A-to-I RNA editing in bacteria: a review. Front Microbiol. 2023;14 :1143929. doi: 10.3389/fmicb.2023.1143929 36960293
[7] Duan Y, Li H, Cai W. Adaptation of A-to-I RNA editing in bacteria, fungi, and animals. Front Microbiol. 2023;14 :1204080. doi: 10.3389/fmicb.2023.1204080 37293227
[8] Bazak L, Haviv A, Barak M, et al. A-to-I RNA editing occurs at over a hundred million genomic sites, located in a majority of human genes. Genome Res. 2014;24 (3 ):365–376. doi: 10.1101/gr.164749.113 24347612
[9] Porath HT, Knisbacher BA, Eisenberg E, et al. Massive A-to-I RNA editing is common across the metazoa and correlates with dsRNA abundance. Genome Biol. 2017;18 (1 ):185. doi: 10.1186/s13059-017-1315-y 28969707
[10] Basilio C, Wahba AJ, Lengyel P, et al. Synthetic polynucleotides and the amino acid code. V. Proc Natl Acad Sci USA. 1962;48 (4 ):613–616. doi: 10.1073/pnas.48.4.613 13865603
[11] Eisenberg E, Levanon EY. A-to-I RNA editing - immune protector and transcriptome diversifier. Nat Rev Genet. 2018;19 (8 ):473–490. doi: 10.1038/s41576-018-0006-1 29692414
[12] Ma L, Zheng C, Liu J, et al. Learning from the codon table: convergent recoding provides novel understanding on the evolution of A-to-I RNA editing. J Mol Evol. 2024;92 (4 ):488–504. doi: 10.1007/s00239-024-10190-z 39012510
[13] Alon S, Garrett SC, Levanon EY, et al. The majority of transcripts in the squid nervous system are extensively recoded by A-to-I RNA editing. Elife. 2015;4 :e05198. doi: 10.7554/eLife.05198 25569156
[14] Heraud-Farlow JE, Chalk AM, Linder SE, et al. Protein recoding by ADAR1-mediated RNA editing is not essential for normal development and homeostasis. Genome Biol. 2017;18 (1 ):166. doi: 10.1186/s13059-017-1301-4 28874170
[15] Shoshan Y, Liscovitch-Brauer N, Rosenthal JJC, et al. Adaptive proteome diversification by nonsynonymous A-to-I RNA editing in coleoid cephalopods. Mol Biol Evol. 2021;38 (9 ):3775–3788. doi: 10.1093/molbev/msab154 34022057
[16] Xin K, Zhang Y, Fan L, et al. Experimental evidence for the functional importance and adaptive advantage of A-to-I RNA editing in fungi. Proc Natl Acad Sci USA. 2023;120 (12 ):e2219029120. doi: 10.1073/pnas.2219029120 36917661
[17] Birk MA, Liscovitch-Brauer N, Dominguez MJ, et al. Temperature-dependent RNA editing in octopus extensively recodes the neural proteome. Cell. 2023;186 (12 ):2544–2555 e2513. doi: 10.1016/j.cell.2023.05.004 37295402
[18] Rangan KJ, Reck-Peterson SL. RNA recoding in cephalopods tailors microtubule motor protein function. Cell. 2023;186 (12 ):2531–2543 e2511. doi: 10.1016/j.cell.2023.04.032 37295401
[19] Garrett S, Rosenthal JJ. RNA editing underlies temperature adaptation in K+ channels from polar octopuses. Science. 2012;335 (6070 ):848–851. doi: 10.1126/science.1212795 22223739
[20] Liscovitch-Brauer N, Alon S, Porath HT, et al. Trade-off between transcriptome plasticity and genome evolution in cephalopods. Cell. 2017;169 (2 ):191–202 e111. doi: 10.1016/j.cell.2017.03.025 28388405
[21] Xu G, Zhang J. In search of beneficial coding RNA editing. Mol Biol Evol. 2015;32 (2 ):536–541. doi: 10.1093/molbev/msu314 25392343
[22] Zhao T, Ma L, Xu S, et al. Narrowing down the candidates of beneficial A-to-I RNA editing by comparing the recoding sites with uneditable counterparts. Nucleus (Calcutta). 2024;15 (1 ):2304503. doi: 10.1080/19491034.2024.2304503
[23] Xu Y, Liu J, Zhao T, et al. Identification and interpretation of A-to-I RNA editing events in insect transcriptomes. Int J Mol Sci. 2023;24 (24 ):17126. doi: 10.3390/ijms242417126 38138955
[24] Lin W, Piskol R, Tan MH, et al. Comment on “widespread RNA and DNA sequence differences in the human transcriptome”. Science. 2012;335 (6074 ):1302. doi: 10.1126/science.1210624
[25] Li Q, Wang Z, Lian J, et al. Caste-specific RNA editomes in the leaf-cutting ant Acromyrmex echinatior. Nat Commun. 2014;5 (1 ):4943. doi: 10.1038/ncomms5943 25266559
[26] Li M, Wang L IX, Bruzel Y, et al. Widespread RNA and DNA sequence differences in the human transcriptome. Science. 2011;333 (6038 ):53–58. doi: 10.1126/science.1207018 21596952
[27] Kleinman CL, Majewski J. Comment on “widespread RNA and DNA sequence differences in the human transcriptome”. Science. 2012;335 (6074 ):1302. doi: 10.1126/science.1209658 author reply 1302
[28] Di Giorgio S, Martignano F, Torcia MG, et al. Conticello SG: evidence for host-dependent RNA editing in the transcriptome of SARS-CoV-2. Sci Adv. 2020;6 (25 ):eabb5813. doi: 10.1126/sciadv.abb5813 32596474
[29] Song Y, He X, Yang W, et al. Virus-specific editing identification approach reveals the landscape of A-to-I editing and its impacts on SARS-CoV-2 characteristics and evolution. Nucleic Acids Res. 2022;50 (5 ):2509–2521. doi: 10.1093/nar/gkac120 35234938
[30] Wei L. Reconciling the debate on deamination on viral RNA. J Appl Genet. 2022;63 (3 ):583–585. doi: 10.1007/s13353-022-00698-9 35507138
[31] Rajakumar R, San Mauro D, Dijkstra MB, et al. Ancestral developmental potential facilitates parallel evolution in ants. Science. 2012;335 (6064 ):79–82. doi: 10.1126/science.1211451 22223805
[32] Smith CR, Toth AL, Suarez AV, et al. Genetic and genomic analyses of the division of labour in insect societies. Nat Rev Genet. 2008;9 (10 ):735–748. doi: 10.1038/nrg2429 18802413
[33] Porath HT, Hazan E, Shpigler H, et al. RNA editing is abundant and correlates with task performance in a social bumblebee. Nat Commun. 2019;10 (1 ):1605. doi: 10.1038/s41467-019-09543-w 30962428
[34] Duan Y, Xu Y, Song F, et al. Differential adaptive RNA editing signals between insects and plants revealed by a new measurement termed haplotype diversity. Biol Direct. 2023;18 (1 ):47. doi: 10.1186/s13062-023-00404-7 37592344
[35] Ramaswami G, Zhang R, Piskol R, et al. Identifying RNA editing sites using RNA sequencing data alone. Nat Methods. 2013;10 (2 ):128–132. doi: 10.1038/nmeth.2330 23291724
[36] Sapiro AL, Shmueli A, Henry GL, et al. Illuminating spatial A-to-I RNA editing signatures within the Drosophila brain. P Natl Acad Sci USA. 2019;116 (6 ):2318–2327. doi: 10.1073/pnas.1811768116
[37] Benjamini Y, Hochberg Y. Controlling the false discovery rate - a practical and powerful approach to multiple testing. J Roy Stat Soc B Met. 1995;57 (1 ):289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
[38] Yablonovitch AL, Deng P, Jacobson D, et al. The evolution and adaptation of A-to-I RNA editing. PLOS Genet. 2017;13 (11 ):e1007064. doi: 10.1371/journal.pgen.1007064 29182635
[39] Zheng C, Ma L, Song F, et al. Comparative genomic analyses reveal evidence for adaptive A-to-I RNA editing in insect Adar gene. Epigenetics. 2024;19 (1 ):2333665. doi: 10.1080/15592294.2024.2333665 38525798
[40] UniProt C. UniProt: a hub for protein information. Nucleic Acids Res. 2015;43 (Database issue ):204–212.
[41] Duan Y, Ma L, Liu J, et al. The first A-to-I RNA editome of hemipteran species Coridius chinensis reveals overrepresented recoding and prevalent intron editing in early-diverging insects. Cell Mol Life Sci. 2024;81 (1 ):136. doi: 10.1007/s00018-024-05175-6 38478033
[42] Duan Y, Ma L, Zhao T, et al. Conserved A-to-I RNA editing with non-conserved recoding expands the candidates of functional editing sites. Fly (Austin). 2024;18 (1 ):2367359. doi: 10.1080/19336934.2024.2367359 38889318
[43] Ma L, Zheng C, Xu S, et al. A full repertoire of Hemiptera genomes reveals a multi-step evolutionary trajectory of auto-RNA editing site in insect adar gene. RNA Biol. 2023;20 (1 ):703–714. doi: 10.1080/15476286.2023.2254985 37676051
[44] Savva YA, Jepson JE, Sahin A, et al. Auto-regulatory RNA editing fine-tunes mRNA re-coding and complex behaviour in Drosophila. Nat Commun. 2012;3 (1 ):790. doi: 10.1038/ncomms1789 22531175
[45] Zhang Y, Duan Y. Genome-wide analysis on driver and passenger RNA editing sites suggests an underestimation of adaptive signals in insects. Genes (Basel). 2023;14 (10 ):1951. doi: 10.3390/genes14101951 37895300
[46] Lewontin RC. On measures of gametic disequilibrium. Genetics. 1988;120 (3 ):849–852. doi: 10.1093/genetics/120.3.849 3224810
[47] Adetula AA, Fan X, Zhang Y, et al. Landscape of tissue-specific RNA editome provides insight into co-regulated and altered gene expression in pigs (Sus-scrofa). RNA Biol. 2021;18 (sup1 ):439–450. doi: 10.1080/15476286.2021.1954380 34314293
[48] Porath HT, Schaffer AA, Kaniewska P, et al. A-to-I RNA editing in the earliest-diverging eumetazoan phyla. Mol Biol Evol. 2017;34 (8 ):1890–1901. doi: 10.1093/molbev/msx125 28453786
[49] Yu Y, Zhou H, Kong Y, et al. The landscape of A-to-I RNA editome is shaped by both positive and purifying selection. PLOS Genet. 2016;12 (7 ):e1006191. doi: 10.1371/journal.pgen.1006191 27467689
[50] Xu G, Zhang J. Human coding RNA editing is generally nonadaptive. Proc Natl Acad Sci USA. 2014;111 (10 ):3769–3774. doi: 10.1073/pnas.1321745111 24567376
[51] Zhang F, Lu Y, Yan S, et al. SPRINT: an snp-free toolkit for identifying RNA editing sites. Bioinformatics. 2017;33 (22 ):3538–3548. doi: 10.1093/bioinformatics/btx473 29036410
[52] Zhang J, Xu C. Gene product diversity: adaptive or not? Trends Genet. 2022;38 (11 ):1112–1122. doi: 10.1016/j.tig.2022.05.002 35641344
[53] Wallberg A, Bunikis I, Pettersson OV, et al. A hybrid de novo genome assembly of the honeybee, apis mellifera, with chromosome-length scaffolds. BMC Genomics. 2019;20 (1 ):275. doi: 10.1186/s12864-019-5642-0 30961563
[54] Graveley BR, Brooks AN, Carlson JW, et al. The developmental transcriptome of Drosophila melanogaster. Nature. 2011;471 (7339 ):473–479. doi: 10.1038/nature09715 21179090
[55] Liu J, Zheng C, Duan Y. New comparative genomic evidence supporting the proteomic diversification role of A-to-I RNA editing in insects. Mol Genet Genomics. 2024;299 (1 ):46. doi: 10.1007/s00438-024-02141-6 38642133
[56] Rodriguez J, Menet JS, Rosbash M. Nascent-seq indicates widespread cotranscriptional RNA editing in Drosophila. Mol Cell. 2012;47 (1 ):27–37. doi: 10.1016/j.molcel.2012.05.002 22658416
[57] St Laurent G, Tackett MR, Nechkin S, et al. Genome-wide analysis of A-to-I RNA editing by single-molecule sequencing in Drosophila. Nat Struct Mol Biol. 2013;20 (11 ):1333–1339. doi: 10.1038/nsmb.2675 24077224
[58] Zhang R, Deng P, Jacobson D, et al. Evolutionary analysis reveals regulatory and functional landscape of coding and non-coding RNA editing. PLOS Genet. 2017;13 (2 ):e1006563. doi: 10.1371/journal.pgen.1006563 28166241
[59] Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29 (1 ):15–21. doi: 10.1093/bioinformatics/bts635 23104886
[60] Camacho C, Coulouris G, Avagyan V, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10 (1 ):421. doi: 10.1186/1471-2105-10-421 20003500
[61] Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30 (4 ):772–780. doi: 10.1093/molbev/mst010 23329690
[62] Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25 (14 ):1754–1760. doi: 10.1093/bioinformatics/btp324 19451168
[63] Picardi E, Pesole G. Reditools: high-throughput RNA editing detection made easy. Bioinformatics. 2013;29 (14 ):1813–1814. doi: 10.1093/bioinformatics/btt287 23742983
[64] McKenna A, Hanna M, Banks E, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20 (9 ):1297–1303. doi: 10.1101/gr.107524.110 20644199
[65] Eisenberg E. Bioinformatic approaches for identification of A-to-I editing sites. Curr Top Microbiol Immunol. 2012;353 :145–162.21751095
[66] Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30 (7 ):923–930. doi: 10.1093/bioinformatics/btt656 24227677
[67] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15 (12 ):550. doi: 10.1186/s13059-014-0550-8 25516281
