Research

Dynamic A-to-I RNA editing in response to gut microbiome in honey bees

    • State Key Laboratory of Agricultural and Forestry Biosecurity, MOA Key Laboratory of Pest Monitoring and Green Management, Department of Entomology, College of Plant Protection, China Agricultural University, Beijing 100193, China
    • 1 These authors contributed equally to this work.
Published March 4, 2026. Vol 36 Issue 4, pp. 713-725. https://doi.org/10.1101/gr.280291.124
Download PDF Cite Article Permissions Share
cover of Genome Research Vol 36 Issue 10
Current Issue:

Abstract

ADAR-mediated adenosine-to-inosine (A-to-I) mRNA editing contributes to the proteomic diversity and behavioral complexity of animals. Although recent studies indicate that the gut microbiome influences various aspects of animal behavior, the potential involvement of RNA editing in gut–brain interactions remains unexplored. Here, we perform comparative transcriptomic analyses between heads of germ-free (GF) versus conventional (CV) honey bees. A total of 1528 A-to-I editing sites are identified in honey bee heads, among which nonsynonymous editing sites are overrepresented compared with random expectation. Overall editing levels are significantly downregulated in GF compared with CV, but the Adar gene is not differentially expressed. However, the Adar p.482 Ile > Met auto-recoding site, which is speculated to modulate Adar activity, is identified as a high-confidence differential editing site (DES) with decreased editing level in GF. Quantification of gut microbiota across 12 CV individuals reveals a significant positive correlation between p.482 Ile > Met editing level and Lactobacillus and Bombilactobacillus abundance. Colonization of a single bacterium Lactobacillus or Bombilactobacillus instead of Gilliamella in GF bees successfully restores the Adar p.482 Ile > Met and the global editing level. Our work demonstrates the complex and dynamic transcriptomic diversity exerted by A-to-I RNA editing and discovers the axis of gut–Lactobacillus/Bombilactobacillus–brain–Adar–global RNA editing, providing an exciting scenario that gut microbiomes could impact RNA editing, which might further facilitate phenotypic plasticity of social insects.


Adenosine-to-inosine (A-to-I) mRNA editing, mediated by the deamination enzyme adenosine deaminase acting on RNA (ADAR), is a prevalent type of RNA modification in metazoans (Palladino et al. 2000a; Savva et al. 2012b). Inosine is recognized as guanosine by cellular machineries, and therefore, A-to-I RNA editing leads to similar consequences to A-to-G mutations, diversifying the transcriptome (Eisenberg and Levanon 2018; Duan 2025). Editing events taking place in coding sequences (CDSs) might alter the genomically encoded amino acid and lead to nonsynonymous changes, namely, “recoding” (Alon et al. 2015). Genome-wide profiling of RNA editing sites has been performed in several metazoans, showing us highly distinct distributions of RNA editing events (Liscovitch-Brauer et al. 2017; Porath et al. 2017; Duan et al. 2023a; Zhang et al. 2023). Although extensive RNA editing on noncoding repetitive elements was observed in most species (Porath et al. 2017; Cao and Duan 2025; Duan et al. 2025a), recoding is overrepresented in two major clades, coleoids (Cephalopoda) (Alon et al. 2015; Liscovitch-Brauer et al. 2017) and Drosophila (Insecta) (Yu et al. 2016; Yablonovitch et al. 2017a; Duan et al. 2024a).

The temporospatial regulation of RNA editing circumvents the pleiotropic effects caused by DNA mutations, facilitating the organisms to adapt to changing environments. Experimental evidence in the octopus (Birk et al. 2023; Rangan and Reck-Peterson 2023) and fungi (Xin et al. 2023; Qi et al. 2024) supports the advantage of flexible RNA editing under different biotic or abiotic conditions. Other in silico observational studies revealing such flexibility include the temperature-dependent RNA editing in ectotherms (Rieder et al. 2015; Yablonovitch et al. 2017b; Duan et al. 2024b), tissue- and stage-specific RNA editing in animals (Graveley et al. 2011; Duan et al. 2023c, 2025b; Zheng et al. 2024b), and caste-specific RNA editing in social insects (Li et al. 2014a; Porath et al. 2019). In a broader range of species and conditions, the flexibility of RNA editing has not been systematically investigated.

The western honey bee (Apis mellifera) is a social insect that has acquired the ability to generate phenotypically distinct castes from the same genome. It has been extensively studied for their complex behaviors (Zayed and Robinson 2012). Although it is believed that RNA editing contributes to such phenotypic plasticity in social insects (Li et al. 2014a; Porath et al. 2019), several recent studies have shown that the honey bee gut microbiota can influence the host behavior and gene expression in the brain. For instance, the gut microbes can enhance learning and memory (Zhang et al. 2022a; Zhong et al. 2024) and modify the sugar sensitivity of honey bees (Zheng et al. 2017; Zhang et al. 2022b). Gut microbiota can also affect the amount of host neurotransmitters (Zhang et al. 2022b) and modulate the social network (Liberti et al. 2022) as well as foraging-related behavior in bees (Liberti et al. 2024; Vernier et al. 2024).

Honey bees are a model organism for gut microbiota research (Zheng et al. 2018) and have a simple and host-specific gut microbiota (Guo et al. 2023; Luo et al. 2024). Newly emerged adult worker bees contain almost no gut bacteria and acquire the gut microbes through feeding, contact with feces, and hive materials (Martinson et al. 2012). The gut microbiota of workers in A. mellifera is mainly composed of five core bacterial phylotypes: Gilliamella, Snodgrassella, Lactobacillus Firm-4 (currently reclassified as Bombilactobacillus), Lactobacillus Firm-5 (currently reclassified as Lactobacillus), and Bifidobacterium (Motta and Moran 2024). Other noncore gut-restricted bacteria include Bartonella, Frischella, and Commensalibacter (Motta and Moran 2024). These eight lineages account for >95% of the total gut microbiota content (Kwong and Moran 2016). A phylotype is defined by a 97% similarity in 16S rRNA sequences (Moreira and López-García 2015). These bacterial phylotypes can be cultured in vitro. Under laboratory conditions, germ-free bees can be colonized with gut bacteria to restore the conventional normal state, which provides a good system for studying the functions of gut microbiota such as in the gut–brain axis (Cryan et al. 2019).

A-to-I RNA editing as a means of proteomic diversifier in nervous system affects various neuronal processes and animal behaviors (Li et al. 2014a; Porath et al. 2019), but its connection to gut microbiota remains unclear. In this work, we bridge this gap by comparing the RNA editomes of conventional (CV) and germ-free (GF) honey bees, revealing the complex and dynamic transcriptomic diversity exerted by A-to-I RNA editing, and discover the roles of Adar in the gut microbiome-brain axis. Our study indicates that RNA editing may regulate the honey bee's phenotypic and behavioral complexity through interactions with gut microbiomes.

Results

The new version of A. mellifera genome contains a single Adar protein

The new A. mellifera genome version HAv3.1 was updated in 2019 (Wallberg et al. 2019). Similar to the previous investigation in the bumble bee (Porath et al. 2019), it is essential to confirm the completeness of Adar gene in the new genome. In the A. mellifera genome version HAv3.1, two transcripts (XM_026445265.1 and XM_006563327.3) are annotated. We cloned those transcripts with RACE and only validated one transcript XM_006563327.3 with a shorter 3′ UTR, not affecting the annotated CDS (Supplemental Fig. S1). The full-length cloned Adar transcript is 2545 bp, encoding a predicted protein of 663 amino acids (AAs). The N terminus contains two dsRNA-binding domains (AA positions 78–140 and 197–241), and the C terminus contains one deaminase domain (AA positions 294–659).

A-to-I RNA editomes in CV and GF honey bees

To investigate the effect of the gut microbiome on A-to-I RNA editing profile in A. mellifera, we fed the homogenates of nurse bee hindguts to the newly emerged GF bees to generate CV bees (Methods). Subsequently, we characterized the RNA editomes of CV and GF groups, each comprising three 7-day-old individuals (Fig. 1A). The heads of each group were used to construct an RNA-seq library, and their respective bodies were subjected to DNA resequencing (Fig. 1A). On average, we obtained 11.56 Gb RNA-seq data and 5.78 Gb DNA-seq data for each of the six samples, covering an average genome-wide depth of 26× for DNA-seq and 68× for RNA-seq (for statistics, see Table 1; for expression profile of six samples, see Supplemental Fig. S2). Throughout this paper, CV bees are regarded as the control group, whereas GF bees represent an experimentally modified condition in which newly emerged bees are prevented from receiving homogenates from nurse bees.

Figure 1.

A-to-I RNA editomes of honey bee samples. (A) Scheme of identifying A-to-I RNA editing sites in the honey bee samples. Three conventional (CV) and three germ-free (GF) samples were used. Heads were used to construct RNA-seq libraries, and bodies were subjected to DNA resequencing. The reliable RNA–DNA differences (RDDs) were regarded as RNA editing sites. (B) The number of different types of mismatches in the RDDs. Only the RDDs in the genic region were considered. A-to-G RDDs represent A-to-I RNA editing sites. The numbers of A-to-I editing sites for each functional category are shown in a barplot. (C) Histogram showing the distribution of editing levels of A-to-I RNA editing sites. (D) Comparison of editing levels of nonsynonymous and synonymous sites in CV and GF honey bees. P-values were calculated by Wilcoxon rank-sum tests. (E) GO enrichment of genes with nonsynonymous editing in honey bee samples.

713f01
Table 1.

Sequencing depth of samples used in this study

SampleData size (Gb)Genome-wide coverage (×)
RNA-seqDNA-seqRNA-seqaDNA-seq
CV19.665.675725
CV210.496.456229
CV312.735.737525
GF114.274.968422
GF211.736.026927
GF310.475.876226

[i] aThe RNA-seq coverage was calculated by data size divided by the total length of genic region instead of whole genome size. (CV) conventional bees; (GF) germ-free bees.

Following the methods of our previous work (Duan et al. 2024a), we used meticulous pipelines to identify RNA–DNA differences (RDDs; Methods) (Fig. 1A). A union of 2251 RDDs was identified in genic region. We found that 1528 (68%) of the RDDs belonged to A-to-G variation, representing reliable A-to-I RNA editing sites (Fig. 1B; Supplemental Table S1). The numbers of such editing sites per sample ranged from 1295 to 1474 (for details, see Supplemental Table S2). These A-to-I RNA editing sites identified by comparing the DNA and RNA sequences of matched sample were termed “regular editing sites.” The median editing level of 1528 regular editing sites was 14.3%, and only 8.7% of all sites had an editing level >50% (Fig. 1C). These findings resembled the editing patterns observed in leaf-cutting ants (Hymenoptera), in which the median editing level was 12.6%, and ∼8% of sites had a level exceeding 50% (Li et al. 2014a). Moreover, an average of 399 reads covered each editing site (median = 178), and 56 reads supported the editing event on each site (median = 23), demonstrating sufficient depth for a set of reliable editing sites (Supplemental Fig. S3). The validation of individual RNA editing sites (and the corresponding positions in DNA) was performed by Sanger sequencing in subsequent sections.

In parallel, we also utilized a DNA-free method to identify hyperediting sites in the transcriptome (Porath et al. 2014; Zhang et al. 2017). A total of 5 × 104∼105 A-to-I hyperediting sites were detected in six samples (Supplemental Fig. S4). Notably, the hyperediting methodology is prone to finding clustered editing sites in lowly expressed intergenic or repetitive regions (Porath et al. 2017), with an average read support of about 1.5 for editing events in our data (Supplemental Fig. S4). This result is consistent with the pattern observed in bumble bees, in which many hyperediting sites were identified (Porath et al. 2019).

Functional annotation of 1528 A-to-I editing sites revealed 211 nonsynonymous (Nonsyn) sites, 67 synonymous (Syn) sites, and 1250 intronic sites, whereas no UTR sites were detected (Fig. 1B). Note that exons of genes are generally annotated from the beginning of CDSs, and therefore, UTRs are poorly annotated. The Nonsyn/Syn ratio (211/67 = 3.15) was significantly higher than random expectation (4,492,737/1,986,130 = 2.26 if we change all adenosines to guanosines in the reference genome, P = 2.9 × 10−4, Fisher's exact test). This reflects the signals of positive selection on Nonsyn editing sites. Moreover, the median editing levels were higher for Nonsyn sites than for Syn sites (Fig. 1D). We then examined the Gene Ontology (GO) of genes with Nonsyn editing, finding an enrichment in neuron-related genes like transporters, synapses, and ion channels (Fig. 1E).

Differential editing sites between CV and GF honey bees

To identify bona fide differential RNA editing sites (DESs), three replicates might not be sufficient, and thus, we expanded the data size to increase statistical power. We further obtained a second batch of honey bees and generated head transcriptomes of 12 CV and 12 GF samples (Methods). The RNA editing sites identified from the first batch of data were used as anchors. RNA editing levels were quantified for each of the new CV and GF sample or pooled samples. We first compared the global editing levels between CV and GF. We found that CV had significantly higher editing levels compared with GF (Fig. 2A).

Figure 2.

Differential RNA editing between CV and GF. (A) Boxplots showing the overall distribution of pooled editing levels in CV and GF samples. The P-value was calculated by paired Wilcoxon rank-sum test. (B) Volcano plots for differential expression analysis on 12 CV and 12 GF samples. The x-axis represents the log2 fold-change (GF/CV) in gene expression. The y-axis represents the log10-transformed adjusted P-value. The orange dashed line represents the threshold of significance. The Adar gene is labeled in red. (C) Barplots showing Adar expression levels measured by qPCR in 12 CV and 12 GF samples. The P-value was calculated by a t-test. (D) Procedures of defining differential editing sites (DESs). A Wilcoxon rank-sum test comparing 12 CV and 12 GF samples was applied to the editing level of each site. Sites with P < 0.05 were defined as DESs. Fisher's exact tests were then applied to the pooled allele counts of CV and GF. DESs with P < 0.05 were further defined as stringent DESs. (E) Scatterplot showing the editing levels in CV and GF. Levels from pooled reads were used. Colored sites are DESs belonging to different genomic categories. (F) Pie charts showing the annotation of DESs and heatmap illustrating the editing levels of stringent DESs with CV > GF. (G) Functional enrichment of significantly up- and downregulated genes in GF. Circle size represents number of genes. (BP) Biological process; (MF) molecular function; (CC) cellular component.

713f02

Adar is the sole factor that mediated A-to-I RNA editing in insects (Palladino et al. 2000a). The overall difference in editing efficiency is attribute to either the differential expression of Adar gene or the altered catalytic activity of Adar enzyme. To investigate these two possibilities, we first performed differential expression analysis between CV and GF. Totally 1144 differentially expressed genes (DEGs) were obtained, including 654 significantly upregulated and 490 downregulated genes in GF (FDR < 0.05 and log2|fold-change| > 1; Methods). However, Adar expression did not show significant differences between CV and GF samples (Fig. 2B), and this insignificance was confirmed by qPCR on Adar expression (Fig. 2C). It suggests that the global differential editing level might be mediated by the change in Adar activity, which will be discussed in subsequent sections.

Next, we identified individual RNA editing sites that show remarkable editing level changes between GF and CV. Statistical significance was assessed using comparison of editing levels in the 12 CV and 12 GF samples (Methods) (Fig. 2D). We totally obtained 145 DESs between the two conditions. Because of the overall lower editing levels in GF, 122 (84.1%) of DESs showed downregulation in GF compared with CV, and only 23 sites (15.9%) of DESs were upregulated in GF (Fig. 2E). The 122 downregulated sites included 40 Nonsyn, 10 Syn, and 72 intronic sites, whereas the 23 upregulated sites included six Nonsyn, two Syn, and 15 intronic sites (Fig. 2F). Considering that the comparison of editing levels between the 12 CV and 12 GF samples might identify some sites with low coverage that do not reflect reliable differences in editing efficiency, we further refined the DESs by requiring sufficient sequencing depth in the pooled samples (P < 0.05 under Fisher's exact tests) (Fig. 2D). The stringent criteria yielded 59 DESs, including 49 downregulated (83.1%) and 10 upregulated sites in the GF samples. Among the downregulated sites under GF conditions, 19 were Nonsyn, four were Syn, and 26 were located in introns (Fig. 2F).

We wonder whether differential editing is correlated with the differential expression of host genes. We compared the expression level (reads per kilobase per million mapped reads [RPKM]) of the edited genes in CV and GF. The 140 downregulated sites in GF mapped to 81 unique genes, and these genes did not show differential RPKM values between the two conditions (P = 0.12, Wilcoxon rank-sum test). Similarly, the 15 upregulated sites in GF mapped to 12 unique genes, with no significant difference in expression level between two conditions (P = 0.48, Wilcoxon rank-sum test). This suggests that the differential RNA editing levels at specific sites are more likely to reflect changes in editing activity rather than being a consequence of altered host gene expression.

Moreover, because the removal of RNA editing will trigger autoimmune response (Liddicoat et al. 2015), we question whether there are inflammatory signatures in GF bees in which RNA editing is overall decreased. We performed functional enrichment of the 654 significantly upregulated and 490 downregulated genes in GF bees, finding that the upregulated genes indeed showed autophagy signals (Fig. 2G). This observation aligns well with the lower editing levels in GF.

Variation in microbial composition in CV bees correlates with the greater molecular variations in the host compared to GF bees

We sought to validate the DESs using Sanger sequencing. Prior to this, we observed a general trend in which CV samples exhibited greater molecular variation compared with GF samples. Principal component analysis (PCA) based on both gene expression and RNA editing levels revealed that CV samples were widely dispersed, whereas GF samples formed a tightly clustered group (Fig. 3A), suggesting greater intragroup variation among CV samples. Consistently, the standard deviation of RNA editing levels was significantly greater in CV samples than in GF samples (Fig. 3B).

Figure 3.

Intragroup variations of gut microbiota of CV and GF bees. (A) Principal component analyses of gene expression (RPKM) and editing levels of CV and GF samples. CV samples are in orange, and GF samples are in green. (B) Boxplot comparing the standard deviation of editing level in CV and GF. P-value was calculated by Wilcoxon rank-sum test. (C) Distribution of total microbiota abundance of 12 CV and 12 GF samples. The y-axis is log10-transformed. (D) Microbial lineage composition in each of the 12 CV samples. (E) Barplots showing the numbers of editing sites and genes with significant correlation (adjusted P < 0.1) with the relative abundance of each microbial lineage. (F,G) Spearman's correlation between the editing level of the Adar p.482 Ile > Met site and the relative proportion of Lactobacillus (F) and Bombilactobacillus (G) in each CV bee.

713f03

This phenomenon could potentially be explained by variations in gut microbiota composition among the CV bees. To test this hypothesis, we quantified the total microbial abundance across all samples and found that GF bees harbored negligible microbiota (traditionally defined as 105 CFUs), whereas CV bees exhibited high and variable microbial loads, ranging from 109 to 1010 CFUs (Fig. 3C). 16s rRNA sequencing data from bee gut further revealed substantial variation in the relative abundance of microbial lineages across the 12 CV samples (Fig. 3D). The most abundant microbe genera are Bartonella, Snodgrassella, and Lactobacillus.

We hypothesized that the variations in gene expression or RNA editing levels in the CV samples might be attributed to specific microbial lineage. We performed a correlation analysis between the proportion of each microbe in individual CV samples and the corresponding expression or editing levels. In total, 40 unique microbial lineages were identified, among which 15 microbes were detected in at least four out of the 12 CV individuals. We found that Lactobacillus and Bartonella had the largest number of correlated editing sites, whereas Sphingomonas had the largest number of correlated genes (for results of adjusted P < 0.1, see Fig. 3E; for results of raw P < 0.05, see Supplemental Fig. S5). As expected, Adar expression showed no significant correlation with any microbial lineage (adjusted P-values ranged from 0.79 to one for different microbes), further supporting the notion that microbiota changes do not drive differential Adar expression. Next, we examined the overlap between stringent DESs and those sites showing significant correlation with Lactobacillus and Bartonella, the genera selected based on their high abundance (Fig. 3D) and remarkable correlation with editing sites among all microbial lineages (Fig. 3E). It turns out that only one site was obtained, the Adar p.482 Ile > Met auto-recoding site. The editing levels at this site in 12 CV bees were positively correlated with the relative abundance of Lactobacillus (Fig. 3F) and Bombilactobacillus (Fig. 3G) in each individual. Notably, both genera were the core bacteria in honey bees. Bombilactobacillus was historically classified within Lactobacillus and has only recently been redefined as a separate genus (Zheng et al. 2020; Luo et al. 2024), and their close phylogenetic relationship suggests substantial functional similarity between the two groups.

Sanger validation of the differential editing of the Adar p.482 Ile > Met site

Insect Adar protein has two dsRNA-binding domain and one catalytic deaminase domain, and the p.482 Ile > Met auto-recoding site is located in deaminase domain (Fig. 4A; Duan et al. 2023b). Previous studies have reported that the p.482 Ile > Met editing site is broadly conserved across multiple Hymenoptera species (Duan et al. 2023b; Zheng et al. 2024a) and that its auto-editing level correlates with global editing activity in bumble bees (Porath et al. 2019), suggesting a regulatory role of Ile > Met editing in modulating Adar activity. Consistent with this, our results showed that Adar p.482 Ile > Met editing levels were significantly correlated with per-sample cumulative editing levels (a proxy for global editing activity) across CV and GF individuals (Fig. 4B). Moreover, CV samples exhibited significantly higher editing levels at this auto-recoding site compared with GF samples as revealed by our DES analysis (Fig. 4C). The differential editing at the Adar p.482 Ile > Met site was further confirmed by manual inspection and visualization of transcriptome alignment (Fig. 4D) and was validated by Sanger sequencing of the 12 CV and 12 GF samples (Fig. 4E,F).

Figure 4.

Demonstration and validation of differential editing at Adar p.482 Ile > Met site between CV and GF. (A) Protein architecture of Adar. The p.482 Ile > Met recoding site is highlighted. (B) Pearson's correlation between the Adar p.482 Ile > Met editing level and the overall cumulative editing level in each sample. (C) Boxplot comparing the Adar p.482 Ile > Met editing levels in 12 CV and 12 GF transcriptomes (NGS data). The P-value was calculated by Wilcoxon rank-sum test. (D) Integrative Genomics Viewer (IGV) (Robinson et al. 2011) visualization of the Adar p.482 Ile > Met editing site in DNA-seq and RNA-seq. Genome resequencing was shown for the first batch data (three CV and three GF), and the transcriptome was shown for the second batch data without matched DNA-seq (12 CV and 12 GF). A randomly chosen sample for each group was displayed. (E) Sanger sequencing of the Adar p.482 Ile > Met site. Sanger traces of representative DNA and RNA samples were shown. (F) Boxplot showing the Sanger editing levels of the Adar p.482 Ile > Met site in 12 CV and 12 GF samples. The P-value was calculated by Wilcoxon rank-sum test.

713f04

This robust pattern raises the possibility that the core gut bacteria Lactobacillus/Bombilactobacillus might regulate the Adar p.482 Ile > Met auto-editing level in the head, thereby affecting the Adar editing activity and ultimately shaping the global neuronal editome.

Colonization of Lactobacillus/Bombilactobacillus in GF restores the Adar Ile > Met editing level as well as global editome

To test this Lactobacillus/Bombilactobacillus–Adar–editome axis, we downloaded the brain transcriptomes generated by previous studies that align with our research requirements (Zhang et al. 2022a,b). These data include the brain transcriptomes from CV and GF bees, together with the GF bees inoculated with a single bacterial species from Gilliamella (Gi), Bombilactobacillus (Firm-4 [F4]), or Lactobacillus (Firm-5 [F5]), and three biological replicates were performed for each condition (Fig. 5A). According to the expressional analyses of the original literatures, Adar was not differentially expressed between GF and any other conditions (Zhang et al. 2022a,b). If Lactobacillus and Bombilactobacillus are the primary bacterial groups affecting the Adar p.482 Ile > Met and global editing levels, we should expect a differential editing pattern between GF and any of the CV, F4, and F5 samples, whereas GF and Gi should exhibit no substantial difference.

Figure 5.

Validation of the gut bacteria–Adar–editome axis using brain transcriptomes of single bacterial inoculation. (A) Experimental design of the public data. The CV and GF bees were defined similarly to our study. The monocolonized bees were generated by inoculating a single microbial species (F4, F5, or Gi) to the GF bees, respectively. (B) Comparison of global RNA editing levels. Reads from three biological replicates were pooled. Wilcoxon rank-sum tests were used to determine the significance. (C) Comparison of the Adar p.482 Ile > Met editing level. Reads from three biological replicates were pooled. Fisher's exact tests were used to determine the significance. (D) IGV visualization of the Adar p.482 Ile > Met editing level in public transcriptomes. (E) Venn diagram showing the overlap of downregulated nonsynonymous editing sites in GF from different comparisons. (*) P < 0.05; (**) P < 0.01; (***) P < 0.001; (n.s.) not significant.

713f05

We quantified the editing status on 1528 candidate editing sites, finding that the overall editing levels were higher for the CV, F4, and F5 samples, whereas Gi and GF samples were lower (Fig. 5B). Particularly, we focused on the Adar p.482 Ile > Met site and observed that the editing level was higher in CV, F4, and F5 and lower in Gi and GF (Fig. 5C). This trend was confirmed by manual inspection of the sequence alignment (Fig. 5D).

Differential editing analysis revealed 208 DESs between CV and GF (39 up and 169 down in GF), 255 DESs between F4 and GF (81 up and 174 down in GF), and 244 DESs between F5 and GF (78 up and 166 down in GF). We then compared the nonsynonymous sites downregulated in GF across these data sets with the DESs from our 12 CV versus 12 GF samples (Fig. 5E). Four nonsynonymous sites were shared, but the Adar p.482 Ile > Met site was not included. The exclusion resulted from the marginally significant P-value for this site in the public data set of (three CV vs. three GF) under a two-tailed Fisher's exact test (Fig. 5C). Notably, when one-tailed test was applied, the Adar p.482 Ile > Met site showed significant difference across all comparisons.

Validation of Nonsyn editing sites by Sanger sequencing and the functional implication of their differential editing

To better understand the specific functional divergence of individual DES under different conditions, we selected the highly reliable DESs and performed Sanger validation. The Adar p.482 Ile > Met site has been validated in Figure 4. We focused on the four significantly downregulated nonsynonymous sites in GF bees shared by different comparisons (Fig. 5E).

Honey bee transcript XM_006563827.3 encodes a phosphatidylinositol-binding clathrin assembly protein (PICALM) LAP. PICALM is typically involved in clathrin-mediated endocytosis of protein trafficking, and the mutations in this gene are associated with Alzheimer's disease (AD) and other neuropathological phenotypes (Schubert et al. 2012; Merthan et al. 2019), implying the gene's crucial function in maintaining neuronal homeostasis. This LAP gene has a high-confidence DES p.464 Arg > Gly downregulated in GF. The differential editing is validated by Sanger sequencing (Fig. 6A).

Figure 6.

Sanger validation of high-confidence nonsynonymous DESs with significantly lower editing levels in GF compared to other conditions. Sanger traces for representative samples were displayed. Boxplots show the comparison of Sanger editing levels between 12 CV and 12 GF. (A) Phosphatidylinositol-binding clathrin assembly protein (PICALM) LAP. (B) Split ends protein spen. (C) Ly-6 protein qvr. Note that 12 CV and 12 GF samples were subjected to PCR, but only nine CV and six GF samples were successfully sequenced. (D) Voltage-gated potassium channel protein Shab and leukocyte receptor member 8.

713f06

Transcript XM_026439421.1 is orthologous to Drosophila gene FBgn0016977 (spen) that encodes split ends protein. Spen is involved in neuronal cell fate and axon guidance during development (Gu et al. 2017). Drosophila lacking spen has deficiency in lipid and protein catabolism, which influences phenotypic responses to diets (Gillette et al. 2020). As validated by Sanger sequencing, honey bee spen has a p.322 Tyr > Cys recoding site with significantly higher editing level in CV than in GF (Fig. 6B).

Transcript XM_006570837.3 is orthologous to Drosophila gene quiver (qvr; also known as sleepless). This gene encodes an extracellular modulator of the potassium channel Shaker (Ruan et al. 2017). Honey bee qvr gene has a p.106 His > Arg recoding site with significantly lower level in GF than all other conditions (Fig. 6C).

Then, voltage-gated potassium channel protein Shab (XM_026443517.1) has a p.62 Thr > Ala recoding site showing a differential editing level (Fig. 6D). Different from the universally conserved Shab Ile > Met recoding site in metazoans (Liscovitch-Brauer et al. 2017; Porath et al. 2019; Duan et al. 2024a), this Thr > Ala recoding appears to be honey bee specific. The different extent of editing in CV and GF might reflect a response to microbiota unique to the honey bee.

In addition to these four nonsynonymous sites shared by different comparisons in which GF bees constantly have significantly lower editing levels (Fig. 5E), we also validated other DESs identified in our 12 CV versus 12 GF comparison. For example, transcript XM_006561721.3 encoding a leukocyte receptor (member 8) has a p.363 Lys > Arg recoding site with a lower level in GF (Fig. 6D). Although these recoding sites are not always located in the preliminarily annotated domains, their regulatory role in protein folding, function, and translocation could not be excluded.

Taken together, our results suggest a model that the presence of gut microbe Lactobacillus/Bombilactobacillus elevates the editing level of the Adar p.482 Ile > Met auto-recoding site, thereby enhancing Adar enzymatic activity and, in turn, modulating the global editome (Fig. 7).

Figure 7.

The Lactobacillus/Bombilactobacillus–Adar–editome axis. The presence of gut microbe Lactobacillus/Bombilactobacillus elevates the editing level of the Adar p.482 Ile > Met auto-recoding site, increasing the enzymatic activity of Adar, and consequently affects the global editome.

713f07

Discussion

The gut–brain axis serves as a two-way communication pathway facilitating interactions between the gut and the brain. Through this axis, the gut microbiota can influence brain functions (Martin et al. 2018; Mayer et al. 2022), playing a regulatory role in neurotransmitter levels, synaptic plasticity, and neurogenesis in the host central nervous system (Cryan and Dinan 2012; Sherwin et al. 2019). Meanwhile, A-to-I mRNA editing mediated by Adar also regulates a wide range of neuron-related traits in animals (Mehler and Mattick 2007; Li et al. 2014b; Rajendren et al. 2021). Disruption of RNA editing will cause multiple neurological disorders in humans (Slotkin and Nishikura 2013; Yamashita and Kwak 2014; Song et al. 2022) and other mammals (Sommer et al. 1991; Egebjerg and Heinemann 1993; Seeburg et al. 1998; Higuchi et al. 2000), as well as abnormal behaviors in insects and other invertebrates (Palladino et al. 2000b; Tonkin et al. 2002; Savva et al. 2012a). The contribution of RNA editing to complex behavior is especially noteworthy in eusocial insects owing to their ability to generate phenotypically diversified individuals under the same genome (Li et al. 2014a; Porath et al. 2019).

In this work, we investigated the head/brain transcriptomes of CV and GF honey bees. We found that the Adar p.482 Ile > Met auto-recoding level, together with the overall editing levels, was downregulated in GF. However, because Adar itself was not differentially expressed between CV and GF, it raises a compelling question of how the global RNA editing profile is differentially regulated between the two conditions. It is previously speculated that this Ile > Met editing site in the bumble bee positively regulates the enzymatic activity of Adar (Porath et al. 2019). The observation that the Ile > Met recoding level is positively correlated with the transcriptome-wide editing level across different individuals supports this assumption.

Although the gut bacterial composition of adult worker honey bees is relatively stable compared with that of other insects (Engel and Moran 2013), notable individual variations exist among members of the same colony (Moran et al. 2012; Rothman et al. 2018). The honey bee gut microbiota exhibits spatiotemporal dynamics and is influenced by host genetics, nutrition, pesticide exposure, and pathogens (Luo et al. 2024). In the second batch of samples collected from Sanya, Hainan, China, we observed a notably high abundance of Bartonella in many individuals, which may be associated with the sampling season and location. Among the differential RNA editing sites identified in the 12 GF and 12 CV bees, a moderate number overlapped with those identified from public data sets of GF and CV samples. This slight discrepancy may be attributed to variations in the gut microbial composition among different CV bees. Although microbial colonization consistently increased RNA editing levels at the Adar p.482 Ile > Met site and elevated global editing activity, only a subset of editing sites remained stable across conditions. The majority of editing events appeared to vary with changes in microbiota composition. These RNA editing sites are not entirely fixed but exhibit a certain degree of context dependency and plasticity.

Our study suggests the distinct roles for different gut microbes in the gut–brain–RNA editing axis and provide further evidence supporting a potential regulatory role of Lactobacillus and Bombilactobacillus in Adar auto-editing, which in turn affects the global RNA editome. The inoculation of these bacterial species might alter the neurotransmitter level, sucrose responsiveness, and learning efficiency in honey bees (Zhang et al. 2022a,b).

The precise mechanisms through which gut microbes affect Adar auto-editing also need further elucidation. Current evidence suggests that small-molecule metabolites derived from gut microbiota can enter the circulatory system and participate in gut–brain communications (Krautkramer et al. 2021; Ahmed et al. 2022). These metabolites are likely among the key factors influencing Adar-mediated RNA editing, particularly at the p.482 Ile > Met site. Future studies could compare metabolomes of CV and GF bee groups to identify key differential bacterial metabolites, thereby clarifying the regulatory pathway from specific gut microbes to Adar p.482 Ile > Met RNA editing. Promising candidate metabolites include tryptophan-derived compounds such as indole derivatives, which are produced by Lactobacillus in honey bees and influence learning and memory (Zhang et al. 2022a). Short-chain fatty acids (SCFAs) produced by the bee gut microbiota are also candidate metabolites (Zheng et al. 2017), and they can modulate neurogenesis and behavior in mammals (Koh et al. 2016).

Thus, A-to-I RNA editing in animals is more flexible than we previously understood. Eusocial insects are typical examples showing the powerful epigenetic plasticity that can generate phenotypically distinct individuals from the same genome. RNA editing is such an epigenetic approach that shapes the phenotypic plasticity and animals, especially the eusocial insects (Li et al. 2014a; Porath et al. 2019). In our current work, we further propose a new layer of plasticity conferred by RNA editing that can also respond to the gut microbiome. It suggests that RNA editing is widely used by the eusocial insects to cope with external challenges and meanwhile meet the demand for achieving sociality. Moreover, in future works, it will be interesting to investigate the role of RNA editing in the gut–brain axis among more social insects or other animals. Do other target editing sites respond to the gut microbiome? Does the gut microbiome manipulate the Adar pathway to broadly regulate insect behavior? What benefit do microbes derive from their influence over the RNA editing pathway? Answers to these questions will help us establish a more comprehensive functional landscape of A-to-I RNA editing in the light of adaptive evolution and complete the complex connectome of gut–brain interaction.

Taken together, using honey bees as a model, we delicately demonstrate the complex and dynamic transcriptomic diversity conferred by A-to-I RNA editing. We bridge the gut–brain axis with RNA editing, providing innovative perspectives on how the gut microbiome influences gene regulation in the nervous system.

Methods

Generation of GF and CV honey bees

The honey bee bacterial colonization process followed the protocol of Powell et al. (2014). Late-stage pupae were taken from the brood frames, placed in sterile plastic bins, and then kept in an incubation chamber at 35°C, and these newly emerged bees were germ-free (Kwong and Moran 2016). These GF bees aged 0–24 h were starved for 2 h before being inoculated. The bacterial solution was prepared from the hindguts of 20 nurse bees by homogenization using a tissue grinder (Tiangen OSE-Y50), followed by suspension in 3 mL of PBS.

The first batch of honey bees was collected from an apiary located in Beijing, China, in July 2019. The CV groups were given a mixture of the bacterial solution and an equal volume of sterile sucrose syrup (50%, w/v), whereas the control GF group was fed a mixture of PBS and an equal volume of 50% sterile sucrose syrup. Both groups were also provided with sterilized pollen ad libitum. After 24 h, all bees were fed with 50% sucrose and sterile pollen and kept at 35°C. At 7 days after bacterial colonization in CV bees, the heads and thoraxes of the honey bees were collected and stored at –80°C.

To investigate the relationship between gut bacterial composition and RNA editing, we collected a second batch of honey bees from an apiary in Sanya, Hainan Province, China, in March 2025. To ensure that each honey bee in the CV group was successfully inoculated, we hand-fed each bee 10 µL of a mixture containing bacterial solution and an equal volume of sterile sucrose syrup (50% w/v). The GF control group was hand-fed 10 µL of a mixture of PBS and an equal volume of sterile sucrose syrup (50% w/v). Both groups were provided with sterilized pollen ad libitum. After 24 h, all bees were fed 50% sucrose syrup and sterile pollen and maintained for 7 days at 35°C.

Extraction and sequencing of DNA and RNA from GF and CV honey bees

For the first batch of honey bees, RNA was extracted from the heads of three honey bees using TRIzol and designated as one biological replicate. DNA from the thoraxes of the corresponding honey bees in each biological replicate was extracted. For the second batch of 12 GF and CV bees, the head of each bee was set as one biological replicate and used for RNA sequencing and quantitative PCR (qPCR). The gut of each honey bee was collected for 16S rRNA sequencing. The detailed experimental procedures regarding the first and second batches of sequencing data are described in the Supplemental Methods.

Acquisition of public data

The reference genome and annotation files of honey bee A. mellifera were downloaded from NCBI (https://www.ncbi.nlm.nih.gov/), genome version HAv3.1, submitted by Wallberg et al. (2019). The public brain transcriptomes of CV, GF, F4, F5, and Gi bees were downloaded from the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession numbers PRJNA783076 and PRJNA743412.

Mapping and variant calling

We mapped the DNA-seq reads to the reference genome with BWA-MEM version 0.7.17-r1198 (Li and Durbin 2009). The PCR duplicates were removed using Picard version 2.27.5. For RNA-seq reads, we used two aligners to map the reads to the reference genome. Apart from BWA, we used the two-pass mode of STAR version 2.7.6a (Dobin et al. 2013) to map the transcriptome reads. The PCR duplicates were removed as well. For RNA-seq data, we additionally applied the SplitNCigarReads pipeline in GATK version 4.3.0.0 (McKenna et al. 2010) to facilitate variant calling near splicing junctions. Then, for both DNA-seq and RNA-seq, variants were called against reference genome using GATK HaplotypeCaller version 4.3.0.0 (McKenna et al. 2010). To minimize the confounding effect of sequencing errors (which likely occur at the ends of reads), only bases with a quality greater than 30 were used, which guarantees a sequencing accuracy of at least 99.9% (Xu et al. 2023). We only kept the RNA-seq variants supported by both aligners (BWA and STAR) to reduce misalignments.

Variant annotation

Variants were recorded in variant calling format (VCF) and subsequently annotated with SnpEff version 4.5 (Cingolani et al. 2012). VCF contains the location (chromosome and position) and mutation (reference allele/reads count and alternative allele/reads count) information of each variant. To annotate each variant, the reference genome sequence and annotation of A. mellifera were used build an index (SnpEff build) for this species. Then, SnpEff will assign each variant to a genomic region (e.g., intergenic, exonic, or intronic) according to the genome annotation and the location of the variant. For CDS variants, information is provided on whether this mutation leads to synonymous or nonsynonymous changes, as well as the specific codons and amino acids before and after the mutation.

Identification of A-to-I RNA editing sites

The detailed bioinformatic procedures and parameters of identification of RNA editing sites are described in the Supplemental Methods. In brief, de novo identification of RNA editing sites requires more stringent criteria compared with the quantification of RNA editing status on known editing sites (Xu et al. 2023). The de novo identification method considers the sequencing coverage and number of alternative reads counts of both RNA-seq and matched DNA-seq from the same individual, whereas the quantification of RNA editing on known sites only requires the “pileup” of RNA-seq reads (Supplemental Methods).

Identification of DESs

DESs are referred to the RNA editing sites showing significantly different editing levels between the CV and GF samples. First, Fisher's exact tests were applied to the reference reads count and alternative reads count of pooled CV and pooled GF. If a P-value < 0.05, it suggests a significantly different proportion of reference and alternative reads count between CV and GF. These sites were regarded as DESs.

The pooled-read strategy might ignore the within-group variation, and this variation is especially obvious at the sites with low sequencing coverage. Therefore, we further performed Wilcoxon rank-sum tests to compare the 12 CV and 12 GF individuals, requiring a P-value < 0.05. Editing sites passing both the pooled-read comparison and cross-individual comparison were regarded as stringent DESs.

Nevertheless, the extent of stringency of our criteria does not affect our final conclusion that the global RNA editing activity is reduced in GF compared with CV. Moreover, the editing levels of the specific editing sites we are interested in were verified by Sanger sequencing, and the level alterations between two conditions were manually scrutinized. These key results were not affected by the global pipeline for detecting DESs.

Gene expression analysis in RNA-seq data

The software featureCounts (v 2.0.3) was used to quantify the reads mapped to each gene (Liao et al. 2014). Uniquely mapped reads were used. Differential expression analysis was performed using DESeq2 v1.34.0 (Love et al. 2014). PCA was done by R function “princomp.” Gene expression measured by RPKM was used as input for PCA.

16S rRNA sequence analysis

After demultiplexing, the 16s rRNA sequences were filtered using fastp (v0.19.6) (Chen et al. 2018) and merged using FLASH (v1.2.11) (Magoč and Salzberg 2011). Reads with a quality score below 20, incorrect primers, or ambiguous bases were discarded. High-quality sequences were then denoised using the DADA2 plugin (Callahan et al. 2016) within the QIIME2 pipeline with default parameters (Bolyen et al. 2019). Denoised sequences generated by DADA2 are referred to as amplicon sequence variants (ASVs). To reduce the impact of sequencing depth on downstream analyses, the number of sequences per sample was rarefied to 39,712, which still provided an average Good's coverage of 99.99%. Taxonomic assignment of ASVs was performed using the naïve Bayesian classifier (Wang et al. 2007) integrated in QIIME2, with the SILVA 16S rRNA database (v138.2) as the reference (Quast et al. 2013).

Validation of RNA editing sites by Sanger sequencing

For the p.482 Ile > Met editing site in Adar, we conducted amplification of Adar from both gDNA and cDNA extracted from the same sample (the first batch data of three CV and three GF). For cDNA synthesis, 1 µg of total RNA was incubated with gDNA wiper (Vazyme Biotech) for 2 min at 42°C to remove genomic DNA prior to the reverse-transcription reaction. The cDNA was synthesized using the HiScript III RT supermix (Vazyme Biotech), with the reaction carried out for 15 min at 37°C followed by 5 sec at 85°C. For other editing sites to be validated, the procedures were the same except that cDNA was synthesized with HiScript IV all-in-one ultra RT supermix for qPCR (Vazyme Biotech R433-01) in 12 CV and 12 GF. Sanger sequencing files are provided in Supplemental Data S1.

The 50 µL PCR reaction mixture included 1 µL of cDNA or gDNA template, 20 µM of each primer (see Supplemental Table S3), and 2 × phanta max master mix (dye plus) from Vazyme Biotech. The amplification conditions were denaturing for 30 sec at 94°C; then 30 cycles of denaturing for 15 sec at 94°C, annealing for 15 sec at 58°C, and extension for 15 sec at 72°C; and final extension for 5 min at 72°C. Analysis of the sequencing chromatograms of the PCR products derived from DNA templates confirmed the absence of SNP at the site. The PCR products obtained from cDNA templates were utilized to determine the RNA editing levels of G/(A + G) by assessing the peak heights of base A/G at the editing site.

qPCR analysis

The cDNA for qPCR analysis of Adar was synthesized from 1 µg of total RNA using HiScript III RT supermix (Vazyme Biotech) in a 20 µL reaction volume and subsequently diluted fivefold. Each 10 µL PCR reaction mixture contained 1 µL of cDNA template, 2 µM of each primer (refer to Supplemental Table S3), and 5 µL of Taq Pro universal SYBR qPCR master mix from Vazyme Biotech. The thermal cycling conditions were as follows: initial denaturation for 30 sec at 95°C, followed by 40 cycles of denaturation for 10 sec at 95°C, and annealing/extension for 30 sec at 60°C. Melting curves were generated by heating for 15 sec to 95°C, cooling for 60 sec to 60°C, and then reheating for 15 sec to 95°C. The honey bee actin gene was used as the reference gene with primers sourced from Zufelato et al. (2004), and relative expression levels were determined using the 2–ΔΔCT method (Schmittgen and Livak 2008). Each sample was run in triplicate on the same plate, with a negative control using RNase-free water as the template.

RACE of Adar gene in honey bee

The full-length nucleotide sequence of Adar was amplified using the rapid-amplification of cDNA ends (RACE) technique with the HiScript-TS 5′/3′ RACE kit (Vazyme RA101). Gene-specific primers for 5′ RACE (5′ GSP, 5′ nested GSP) and 3′ RACE (3′ GSP, 3′ nested GSP) were designed (Supplemental Table S3) based on the known sequences of Adar (gene ID: 408529) in NCBI. The detailed experimental procedures are described in the Supplemental Methods.

Protein architecture analysis

Domain annotation of proteins was predicted by InterProscan v5.0 with default parameters (Jones et al. 2014). The protein sequences were used as input files.

Statistical tests

Statistical tests and graphical works were performed in R studio (R version 3.6.3) (R Core Team 2020). For the comparisons involving variables with nonnormal distribution (such as editing levels), nonparameter Wilcoxon rank-sum tests were used.

Data access

The raw sequencing data generated in this study have been submitted to the Genome Sequence Archive (GSA; https://ngdc.cncb.ac.cn/gsa/) under accession numbers CRA018500 (three CV and three GF, head transcriptome, and DNA resequencing of matched individual) and CRA029680 (12 CV and 12 GF head transcriptomes and 12 CV gut microbiota).

Competing interest statement

The authors declare no competing interests.

Acknowledgments

This study is financially supported by the National Natural Science Foundation of China (nos. 32170453 and 32300371), the Ministry of Science and Technology of the People's Republic of China (2023YFD2201804), the Chinese Universities Scientific Fund, the Young Elite Scientist Sponsorship Program by the China Association for Science and Technology (CAST) (no. 2023QNRC001), the Young Elite Scientist Sponsorship Program by the Beijing Association for Science and Technology (BAST) (no. BYESS2023160), and the 2115 Talent Development Program of China Agricultural University. We thank the laboratory members who gave constructive suggestions to this project.

Authors’ contributions: Y. Duan, W.C., X.Z., and S.L. conceptualized and supervised this work. Y. Duan, M.H., L.Y., L.M., Y.W., Y. Du, and S.L. obtained the data used in this study. Y. Duan, Y.W., and J.L. performed data analysis. Y. Duan, M.H., and S.L. wrote the original draft of this manuscript. Y. Duan, M.H., Y.W., Y. Du, J.L., F.S., L.T., W.C., H.L., X.Z., and S.L. revised this manuscript. All authors approved the submission of this manuscript.

Notes

[2] Supplementary material [Supplemental material is available for this article.]

[3] Article published online before print. Article, supplemental material, and publication date are at https://www.genome.org/cgi/doi/10.1101/gr.280291.124.

References

  1. ↵
    Ahmed H, Leyrolle Q, Koistinen V, Kärkkäinen O, Layé S, Delzenne N, Hanhineva K. 2022. Microbiota-derived metabolites as drivers of gut-brain communication. Gut Microbes 14: 2102878. 10.1080/19490976.2022.2102878
  2. ↵
    Alon S, Garrett SC, Levanon EY, Olson S, Graveley BR, Rosenthal JJ, Eisenberg E. 2015. The majority of transcripts in the squid nervous system are extensively recoded by A-to-I RNA editing. eLife 4: e05198. 10.7554/eLife.05198
  3. ↵
    Birk MA, Liscovitch-Brauer N, Dominguez MJ, McNeme S, Yue Y, Hoff JD, Twersky I, Verhey KJ, Sutton RB, Eisenberg E, 2023. Temperature-dependent RNA editing in octopus extensively recodes the neural proteome. Cell 186: 2544–2555.e13. 10.1016/j.cell.2023.05.004
  4. ↵
    Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F, 2019. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol 37: 852–857. 10.1038/s41587-019-0209-9
  5. ↵
    Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. 2016. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13: 581–583. 10.1038/nmeth.3869
  6. ↵
    Cao Q, Duan Y. 2025. RNA modifications as a means of self-recognition and immune protection. Crit Rev Biochem Mol Biol 1–7. 10.1080/10409238.2025.2591027
  7. ↵
    Chen S, Zhou Y, Chen Y, Gu J. 2018. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34: i884–i890. 10.1093/bioinformatics/bty560
  8. ↵
    Cingolani P, Platts A, Wang le L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM. 2012. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6: 80–92. 10.4161/fly.19695
  9. ↵
    Cryan JF, Dinan TG. 2012. Mind-altering microorganisms: the impact of the gut microbiota on brain and behaviour. Nat Rev Neurosci 13: 701–712. 10.1038/nrn3346
  10. ↵
    Cryan JF, O'Riordan KJ, Cowan CSM, Sandhu KV, Bastiaanssen TFS, Boehme M, Codagnone MG, Cussotto S, Fulling C, Golubeva AV, 2019. The microbiota-gut-brain axis. Physiol Rev 99: 1877–2013. 10.1152/physrev.00018.2018
  11. ↵
    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. 2013. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29: 15–21. 10.1093/bioinformatics/bts635
  12. ↵
    Duan Y. 2025. Accurate interpretation of inosines in RNAs: recognized as G is more than basepairing with C. Epigenomics 17: 733–738. 10.1080/17501911.2025.2523237
  13. ↵
    Duan Y, Li H, Cai W. 2023a. Adaptation of A-to-I RNA editing in bacteria, fungi, and animals. Front Microbiol 14: 1204080. 10.3389/fmicb.2023.1204080
  14. ↵
    Duan Y, Ma L, Song F, Tian L, Cai W, Li H. 2023b. Autorecoding A-to-I RNA editing sites in the Adar gene underwent compensatory gains and losses in major insect clades. RNA 29: 1509–1519. 10.1261/rna.079682.123
  15. ↵
    Duan Y, Xu Y, Song F, Tian L, Cai W, Li H. 2023c. Differential adaptive RNA editing signals between insects and plants revealed by a new measurement termed haplotype diversity. Biol Direct 18: 47. 10.1186/s13062-023-00404-7
  16. ↵
    Duan Y, Ma L, Liu J, Liu X, Song F, Tian L, Cai W, Li H. 2024a. 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 81: 136. 10.1007/s00018-024-05175-6
  17. ↵
    Duan Y, Ma L, Zhao T, Liu J, Zheng C, Song F, Tian L, Cai W, Li H. 2024b. Conserved A-to-I RNA editing with non-conserved recoding expands the candidates of functional editing sites. Fly (Austin) 18: 2367359. 10.1080/19336934.2024.2367359
  18. ↵
    Duan Y, Cao Q, Xie Q, Ma L, Cai W, Li H. 2025a. On the origin, evolution, and maintenance of RNA editing. Nucleic Acids Res 53: gkaf947. 10.1093/nar/gkaf947
  19. ↵
    Duan Y, Xie Q, Ma L, Cao Q. 2025b. Abundant A-to-I RNA editing in spermatocytes may suppress transposons to compensate for piRNA downregulation in male germlines. RNA Biol 22: 1–9. 10.1080/15476286.2025.2578573
  20. ↵
    Egebjerg J, Heinemann SF. 1993. Ca2+ permeability of unedited and edited versions of the kainate selective glutamate receptor GluR6. Proc Natl Acad Sci 90: 755–759. 10.1073/pnas.90.2.755
  21. ↵
    Eisenberg E, Levanon EY. 2018. A-to-I RNA editing – immune protector and transcriptome diversifier. Nat Rev Genet 19: 473–490. 10.1038/s41576-018-0006-1
  22. ↵
    Engel P, Moran NA. 2013. The gut microbiota of insects: diversity in structure and function. FEMS Microbiol Rev 37: 699–735. 10.1111/1574-6976.12025
  23. ↵
    Gillette CM, Hazegh KE, Nemkov T, Stefanoni D, D'Alessandro A, Taliaferro JM, Reis T. 2020. Gene-diet interactions: dietary rescue of metabolic effects in spen-depleted Drosophila melanogaster. Genetics 214: 961–975. 10.1534/genetics.119.303015
  24. ↵
    Graveley BR, Brooks AN, Carlson JW, Duff MO, Landolin JM, Yang L, Artieri CG, van Baren MJ, Boley N, Booth BW, 2011. The developmental transcriptome of Drosophila melanogaster. Nature 471: 473–479. 10.1038/nature09715
  25. ↵
    Gu T, Zhao T, Kohli U, Hewes RS. 2017. The large and small SPEN family proteins stimulate axon outgrowth during neurosecretory cell remodeling in Drosophila. Dev Biol 431: 226–238. 10.1016/j.ydbio.2017.09.013
  26. ↵
    Guo L, Tang J, Tang M, Luo S, Zhou X. 2023. Reactive oxygen species are regulated by immune deficiency and Toll pathways in determining the host specificity of honeybee gut bacteria. Proc Natl Acad Sci 120: e2219634120. 10.1073/pnas.2219634120
  27. ↵
    Higuchi M, Stefan M, Single FN, Hartner J, Rozov A, Burnashev N, Feldmeyer D, Sprengel R, Seeburg PH. 2000. Point mutation in an AMPA receptor gene rescues lethality in mice deficient in the RNA-editing enzyme ADAR2. Nature 406: 78–81. 10.1038/35017558
  28. ↵
    Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, McWilliam H, Maslen J, Mitchell A, Nuka G, 2014. InterProScan 5: genome-scale protein function classification. Bioinformatics 30: 1236–1240. 10.1093/bioinformatics/btu031
  29. ↵
    Koh A, De Vadder F, Kovatcheva-Datchary P, Bäckhed F. 2016. From dietary fiber to host physiology: short-chain fatty acids as key bacterial metabolites. Cell 165: 1332–1345. 10.1016/j.cell.2016.05.041
  30. ↵
    Krautkramer KA, Fan J, Bäckhed F. 2021. Gut microbial metabolites as multi-kingdom intermediates. Nat Rev Microbiol 19: 77–94. 10.1038/s41579-020-0438-4
  31. ↵
    Kwong WK, Moran NA. 2016. Gut microbial communities of social bees. Nat Rev Microbiol 14: 374–384. 10.1038/nrmicro.2016.43
  32. ↵
    Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25: 1754–1760. 10.1093/bioinformatics/btp324
  33. ↵
    Li Q, Wang Z, Lian J, Schiøtt M, Jin L, Zhang P, Zhang Y, Nygaard S, Peng Z, Zhou Y, 2014a. Caste-specific RNA editomes in the leaf-cutting ant Acromyrmex echinatior. Nat Commun 5: 4943. 10.1038/ncomms5943
  34. ↵
    Li X, Overton IM, Baines RA, Keegan LP, O'Connell MA. 2014b. The ADAR RNA editing enzyme controls neuronal excitability in Drosophila melanogaster. Nucleic Acids Res 42: 1139–1151. 10.1093/nar/gkt909
  35. ↵
    Liao Y, Smyth GK, Shi W. 2014. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30: 923–930. 10.1093/bioinformatics/btt656
  36. ↵
    Liberti J, Kay T, Quinn A, Kesner L, Frank ET, Cabirol A, Richardson TO, Engel P, Keller L. 2022. The gut microbiota affects the social network of honeybees. Nat Ecol Evol 6: 1471–1479. 10.1038/s41559-022-01840-w
  37. ↵
    Liberti J, Frank ET, Kay T, Kesner L, Monié-Ibanes M, Quinn A, Schmitt T, Keller L, Engel P. 2024. Gut microbiota influences onset of foraging-related behavior but not physiological hallmarks of division of labor in honeybees. mBio 15: e0103424. 10.1128/mbio.01034-24
  38. ↵
    Liddicoat BJ, Piskol R, Chalk AM, Ramaswami G, Higuchi M, Hartner JC, Li JB, Seeburg PH, Walkley CR. 2015. RNA editing by ADAR1 prevents MDA5 sensing of endogenous dsRNA as nonself. Science 349: 1115–1120. 10.1126/science.aac7049
  39. ↵
    Liscovitch-Brauer N, Alon S, Porath HT, Elstein B, Unger R, Ziv T, Admon A, Levanon EY, Rosenthal JJC, Eisenberg E. 2017. Trade-off between transcriptome plasticity and genome evolution in cephalopods. Cell 169: 191–202.e11. 10.1016/j.cell.2017.03.025
  40. ↵
    Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550. 10.1186/s13059-014-0550-8
  41. ↵
    Luo S, Zhang X, Zhou X. 2024. Temporospatial dynamics and host specificity of honeybee gut bacteria. Cell Rep 43: 114408. 10.1016/j.celrep.2024.114408
  42. ↵
    Magoč T, Salzberg SL. 2011. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27: 2957–2963. 10.1093/bioinformatics/btr507
  43. ↵
    Martin CR, Osadchiy V, Kalani A, Mayer EA. 2018. The brain-gut-microbiome axis. Cell Mol Gastroenterol Hepatol 6: 133–148. 10.1016/j.jcmgh.2018.04.003
  44. ↵
    Martinson VG, Moy J, Moran NA. 2012. Establishment of characteristic gut bacteria during development of the honeybee worker. Appl Environ Microbiol 78: 2830–2840. 10.1128/AEM.07810-11
  45. ↵
    Mayer EA, Nance K, Chen S. 2022. The gut-brain axis. Annu Rev Med 73: 439–453. 10.1146/annurev-med-042320-014032
  46. ↵
    McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, 2010. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20: 1297–1303. 10.1101/gr.107524.110
  47. ↵
    Mehler MF, Mattick JS. 2007. Noncoding RNAs and RNA editing in brain development, functional diversification, and neurological disease. Physiol Rev 87: 799–823. 10.1152/physrev.00036.2006
  48. ↵
    Merthan L, Haller A, Thal DR, von Einem B, von Arnim CAF. 2019. The role of PTB domain containing adaptor proteins on PICALM-mediated APP endocytosis and localization. Biochem J 476: 2093–2109. 10.1042/BCJ20180840
  49. ↵
    Moran NA, Hansen AK, Powell JE, Sabree ZL. 2012. Distinctive gut microbiota of honey bees assessed using deep sampling from individual worker bees. PLoS One 7: e36393. 10.1371/journal.pone.0036393
  50. ↵
    Moreira D, López-García P. 2015. Phylotype. In Encyclopedia of astrobiology (ed. Gargaud M, ), p. 1890. Springer, Berlin.
  51. ↵
    Motta EVS, Moran NA. 2024. The honeybee microbiota and its impact on health and disease. Nat Rev Microbiol 22: 122–137. 10.1038/s41579-023-00990-3
  52. ↵
    Palladino MJ, Keegan LP, O'Connell MA, Reenan RA. 2000a. dADAR, a Drosophila double-stranded RNA-specific adenosine deaminase is highly developmentally regulated and is itself a target for RNA editing. RNA 6: 1004–1018. 10.1017/S1355838200000248
  53. ↵
    Palladino MJ, Keegan LP, O'Connell MA, Reenan RA. 2000b. A-to-I pre-mRNA editing in Drosophila is primarily involved in adult nervous system function and integrity. Cell 102: 437–449. 10.1016/S0092-8674(00)00049-0
  54. ↵
    Porath HT, Carmi S, Levanon EY. 2014. A genome-wide map of hyper-edited RNA reveals numerous new sites. Nat Commun 5: 4726. 10.1038/ncomms5726
  55. ↵
    Porath HT, Knisbacher BA, Eisenberg E, Levanon EY. 2017. Massive A-to-I RNA editing is common across the metazoa and correlates with dsRNA abundance. Genome Biol 18: 185. 10.1186/s13059-017-1315-y
  56. ↵
    Porath HT, Hazan E, Shpigler H, Cohen M, Band M, Ben-Shahar Y, Levanon EY, Eisenberg E, Bloch G. 2019. RNA editing is abundant and correlates with task performance in a social bumblebee. Nat Commun 10: 1605. 10.1038/s41467-019-09543-w
  57. ↵
    Powell JE, Martinson VG, Urban-Mead K, Moran NA. 2014. Routes of acquisition of the gut microbiota of the honey bee Apis mellifera. Appl Environ Microbiol 80: 7378–7387. 10.1128/AEM.01861-14
  58. ↵
    Qi Z, Lu P, Long X, Cao X, Wu M, Xin K, Xue T, Gao X, Huang Y, Wang Q, 2024. Adaptive advantages of restorative RNA editing in fungi for resolving survival-reproduction trade-offs. Sci Adv 10: eadk6130. 10.1126/sciadv.adk6130
  59. ↵
    Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO. 2013. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41: D590–D596. 10.1093/nar/gks1219
  60. ↵
    Rajendren S, Dhakal A, Vadlamani P, Townsend J, Deffit SN, Hundley HA. 2021. Profiling neural editomes reveals a molecular mechanism to regulate RNA editing during development. Genome Res 31: 27–39. 10.1101/gr.267575.120
  61. ↵
    Rangan KJ, Reck-Peterson SL. 2023. RNA recoding in cephalopods tailors microtubule motor protein function. Cell 186: 2531–2543.e11. 10.1016/j.cell.2023.04.032
  62. ↵
    R Core Team. 2020. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org/.
  63. ↵
    Rieder LE, Savva YA, Reyna MA, Chang YJ, Dorsky JS, Rezaei A, Reenan RA. 2015. Dynamic response of RNA editing to temperature in Drosophila. BMC Biol 13: 1. 10.1186/s12915-014-0111-3
  64. ↵
    Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. 2011. Integrative genomics viewer. Nat Biotechnol 29: 24–26. 10.1038/nbt.1754
  65. ↵
    Rothman JA, Carroll MJ, Meikle WG, Anderson KE, McFrederick QS. 2018. Longitudinal effects of supplemental forage on the honey bee (Apis mellifera) microbiota and inter- and intra-colony variability. Microb Ecol 76: 814–824. 10.1007/s00248-018-1151-y
  66. ↵
    Ruan H, Ueda A, Xing X, Wan X, Strub B, Mukai S, Certel K, Green D, Belozerov K, Yao WD, 2017. Generation and characterization of new alleles of quiver (qvr) that encodes an extracellular modulator of the Shaker potassium channel. J Neurogenet 31: 325–336. 10.1080/01677063.2017.1393076
  67. ↵
    Savva YA, Jepson JE, Sahin A, Sugden AU, Dorsky JS, Alpert L, Lawrence C, Reenan RA. 2012a. Auto-regulatory RNA editing fine-tunes mRNA re-coding and complex behaviour in Drosophila. Nat Commun 3: 790. 10.1038/ncomms1789
  68. ↵
    Savva YA, Rieder LE, Reenan RA. 2012b. The ADAR protein family. Genome Biol 13: 252. 10.1186/gb-2012-13-12-252
  69. ↵
    Schmittgen TD, Livak KJ. 2008. Analyzing real-time PCR data by the comparative CT method. Nat Protoc 3: 1101–1108. 10.1038/nprot.2008.73
  70. ↵
    Schubert KO, Föcking M, Prehn JH, Cotter DR. 2012. Hypothesis review: are clathrin-mediated endocytosis and clathrin-dependent membrane and protein trafficking core pathophysiological processes in schizophrenia and bipolar disorder? Mol Psychiatry 17: 669–681. 10.1038/mp.2011.123
  71. ↵
    Seeburg PH, Higuchi M, Sprengel R. 1998. RNA editing of brain glutamate receptor channels: mechanism and physiology. Brain Res Brain Res Rev 26: 217–229. 10.1016/S0165-0173(97)00062-3
  72. ↵
    Sherwin E, Bordenstein SR, Quinn JL, Dinan TG, Cryan JF. 2019. Microbiota and the social brain. Science 366: eaar2016. 10.1126/science.aar2016
  73. ↵
    Slotkin W, Nishikura K. 2013. Adenosine-to-inosine RNA editing and human disease. Genome Med 5: 105. 10.1186/gm508
  74. ↵
    Sommer B, Köhler M, Sprengel R, Seeburg PH. 1991. RNA editing in brain controls a determinant of ion flow in glutamate-gated channels. Cell 67: 11–19. 10.1016/0092-8674(91)90568-J
  75. ↵
    Song B, Shiromoto Y, Minakuchi M, Nishikura K. 2022. The role of RNA editing enzyme ADAR1 in human disease. Wiley Interdiscip Rev RNA 13: e1665. 10.1002/wrna.1665
  76. ↵
    Tonkin LA, Saccomanno L, Morse DP, Brodigan T, Krause M, Bass BL. 2002. RNA editing by ADARs is important for normal behavior in Caenorhabditis elegans. EMBO J 21: 6025–6035. 10.1093/emboj/cdf607
  77. ↵
    Vernier CL, Nguyen LA, Gernat T, Ahmed AC, Chen Z, Robinson GE. 2024. Gut microbiota contribute to variations in honey bee foraging intensity. ISME J 18: wrae030. 10.1093/ismejo/wrae030
  78. ↵
    Wallberg A, Bunikis I, Pettersson OV, Mosbech MB, Childers AK, Evans JD, Mikheyev AS, Robertson HM, Robinson GE, Webster MT. 2019. A hybrid de novo genome assembly of the honeybee, Apis mellifera, with chromosome-length scaffolds. BMC Genomics 20: 275. 10.1186/s12864-019-5642-0
  79. ↵
    Wang Q, Garrity GM, Tiedje JM, Cole JR. 2007. Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microbiol 73: 5261–5267. 10.1128/AEM.00062-07
  80. ↵
    Xin K, Zhang Y, Fan L, Qi Z, Feng C, Wang Q, Jiang C, Xu JR, Liu H. 2023. Experimental evidence for the functional importance and adaptive advantage of A-to-I RNA editing in fungi. Proc Natl Acad Sci 120: e2219029120. 10.1073/pnas.2219029120
  81. ↵
    Xu Y, Liu J, Zhao T, Song F, Tian L, Cai W, Li H, Duan Y. 2023. Identification and interpretation of A-to-I RNA editing events in insect transcriptomes. Int J Mol Sci 24: 17126. 10.3390/ijms242417126
  82. ↵
    Yablonovitch AL, Deng P, Jacobson D, Li JB. 2017a. The evolution and adaptation of A-to-I RNA editing. PLoS Genet 13: e1007064. 10.1371/journal.pgen.1007064
  83. ↵
    Yablonovitch AL, Fu J, Li K, Mahato S, Kang L, Rashkovetsky E, Korol AB, Tang H, Michalak P, Zelhof AC, 2017b. Regulation of gene expression and RNA editing in Drosophila adapting to divergent microclimates. Nat Commun 8: 1570. 10.1038/s41467-017-01658-2
  84. ↵
    Yamashita T, Kwak S. 2014. The molecular link between inefficient GluA2 Q/R site-RNA editing and TDP-43 pathology in motor neurons of sporadic amyotrophic lateral sclerosis patients. Brain Res 1584: 28–38. 10.1016/j.brainres.2013.12.011
  85. ↵
    Yu Y, Zhou H, Kong Y, Pan B, Chen L, Wang H, Hao P, Li X. 2016. The landscape of A-to-I RNA editome is shaped by both positive and purifying selection. PLoS Genet 12: e1006191. 10.1371/journal.pgen.1006191
  86. ↵
    Zayed A, Robinson GE. 2012. Understanding the relationship between brain gene expression and social behavior: lessons from the honey bee. Annu Rev Genet 46: 591–615. 10.1146/annurev-genet-110711-155517
  87. ↵
    Zhang F, Lu Y, Yan S, Xing Q, Tian W. 2017. SPRINT: an SNP-free toolkit for identifying RNA editing sites. Bioinformatics 33: 3538–3548. 10.1093/bioinformatics/btx473
  88. ↵
    Zhang Z, Mu X, Cao Q, Shi Y, Hu X, Zheng H. 2022a. Honeybee gut Lactobacillus modulates host learning and memory behaviors via regulating tryptophan metabolism. Nat Commun 13: 2037. 10.1038/s41467-022-29760-0
  89. ↵
    Zhang Z, Mu X, Shi Y, Zheng H. 2022b. Distinct roles of honeybee gut bacteria on host metabolism and neurological processes. Microbiol Spectr 10: e0243821. 10.1128/spectrum.02438-21
  90. ↵
    Zhang P, Zhu Y, Guo Q, Li J, Zhan X, Yu H, Xie N, Tan H, Lundholm N, Garcia-Cuetos L, 2023. On the origin and evolution of RNA editing in metazoans. Cell Rep 42: 112112. 10.1016/j.celrep.2023.112112
  91. ↵
    Zheng H, Powell JE, Steele MI, Dietrich C, Moran NA. 2017. Honeybee gut microbiota promotes host weight gain via bacterial metabolism and hormonal signaling. Proc Natl Acad Sci 114: 4775–4780. 10.1073/pnas.1701819114
  92. ↵
    Zheng H, Steele MI, Leonard SP, Motta EVS, Moran NA. 2018. Honey bees as models for gut microbiota research. Lab Anim (NY) 47: 317–325. 10.1038/s41684-018-0173-x
  93. ↵
    Zheng J, Wittouck S, Salvetti E, Franz C, Harris HMB, Mattarelli P, O'Toole PW, Pot B, Vandamme P, Walter J, 2020. A taxonomic note on the genus Lactobacillus: description of 23 novel genera, emended description of the genus Lactobacillus Beijerinck 1901, and union of Lactobacillaceae and Leuconostocaceae. Int J Syst Evol Microbiol 70: 2782–2858. 10.1099/ijsem.0.004107
  94. ↵
    Zheng C, Liu J, Duan Y. 2024a. Adaptive evolution of A-to-I auto-editing site in Adar of eusocial insects. BMC Genomics 25: 803. 10.1186/s12864-024-10709-0
  95. ↵
    Zheng C, Ma L, Song F, Tian L, Cai W, Li H, Duan Y. 2024b. Comparative genomic analyses reveal evidence for adaptive A-to-I RNA editing in insect Adar gene. Epigenetics 19: 2333665. 10.1080/15592294.2024.2333665
  96. ↵
    Zhong Z, Mu X, Lang H, Wang Y, Jiang Y, Liu Y, Zeng Q, Xia S, Zhang B, Wang Z, 2024. Gut symbiont-derived anandamide promotes reward learning in honeybees by activating the endocannabinoid pathway. Cell Host Microbe 32: 1944–1958.e7. 10.1016/j.chom.2024.09.013
  97. ↵
    Zufelato MS, Lourenço AP, Simões ZL, Jorge JA, Bitondi MM. 2004. Phenoloxidase activity in Apis mellifera honey bee pupae, and ecdysteroid-dependent expression of the prophenoloxidase mRNA. Insect Biochem Mol Biol 34: 1257–1268. 10.1016/j.ibmb.2004.08.005
Loading
Loading
Loading
Loading
Back to top