18 Datasets of IBP infections in mice were analyzed including the lung, liver, right femur, brain, bone-marrow-derived macrophages (BMDM), bone-marrow-derived neutrophils (BMDN), and macrophage cell line raw264.7. All samples from the datasets were firstly combined into an uninfected group and an infected group (Table 1) and subjected to subsequent analysis. Principal component analysis (PCA) showed that the uninfected and infected groups clustered separately based on the editing level of differential RNA editing (DRE) sites (Figure 1A). In terms of the editing level of the DRE sites, most models (13/18) showed higher editing levels after bacterial infection (Figure 1B). The majority of DRE sites in bacterial infections were 3′ -untranslated region (UTR), intronic, and missense variants (Figure 1C). These results suggested distinct alteration of RNA editing profiles during IBP infections.

A-to-I RNA editing is mediated by RNA editing enzymes Adar and Adarb1 (13). Our results showed that Adar expression increased, while Adarb1 decreased in most IBP infections (Figures 2A, B). DRE sites in all infections were further compared to identify shared DRE sites (Supplementary Table S1 and Figure 2C). In particular, DRE sites in Calmodulin 1 (Calm1: chr12: 100207186) and Tyrosine 3-Monooxygenase/Tryptophan 5-Monooxygenase Activation Protein Gamma (Ywhag: chr5: 135909342) were shared by 12 infection datasets, which were also predicted to exert a cis-regulatory effect on the gene expression (Supplementary Figures S2A, B). Gene ontologies (GO) showed that the DRE genes were enriched in immune response pathways, such as neutrophil-mediated immunity and regulation of T cell cytokine production (Supplementary Table S3.), phosphate-containing metabolic compound process, lipid metabolism, and translational regulation (Figure 2D). The KEGG pathway analysis demonstrated enrichment of the DRE genes in lysosomes and autophagy pathways (Supplementary Figure S2C).

Given the high incidence, infectivity, harm (40), and commonality of pneumonia, we further focused on its RNA editing. A-to-I RNA editing was the most frequent among all RNA editing types in terms of both editing sites and edited genes in IBP pneumonia (Supplement Figure S1). Thus we focused on A-to-I RNA editing in subsequent analysis. 656 editing sites in 138 edited genes and 1090 editing sites in 189 edited genes were exclusively detected in uninfected and infected lung tissues of the IBP pneumonia datasets, respectively (Figures 3A, B). Moreover, the number of RNA editing sites and edited genes as well as the editing level showed an up-regulated trend after bacterial infection (Figures 3C–E). More specifically, the top 30 sites that were the most differentially edited in IBP pneumonia datasets were shown in Figure 3F. As shown in Figure 3G, two significant cis-regulatory DRE sites with up-regulated RNA editing levels in infected lung tissues were found in Basic Helix-Loop-Helix Family Member E40 (Bhhe40) (Bhlhe40: chr6:108665779, Spearman r = 0.47, P = 2.7 × 10⁻⁵) and Protein Phosphatase 1 Regulatory Subunit 15B (Ppp1r15b) (Ppp1r15b: chr1:133138010, Spearman r = 0.42, P = 2.6 × 10⁻⁴). Two cis-regulatory DRE sites with down-regulated RNA editing levels were Sideroflexin 2 (Sfxn2) (Sfxn2: chr19:46595684, Spearman r = -0.4, P = 0.009) and Nuclear Factor I A (Nfia) (Nfia: chr4:98118559, with Spearman r = -0.36, P = 0.006) (Figure 3H and Supplementary Table S4). Enrichment analysis and Gene Set Enrichment Analysis (GSEA) revealed that these DRE genes were enriched in functions and pathways related to lipid metabolism, innate immunity, and GTPase-related regulation in IBP pneumonia (see Figures 3I–K).

Comparing the similarities and differences in RNA editing between viral and IBP pneumonia, A-to-I RNA editing was the most frequent among all RNA editing types (Supplementary Figures S1C, D), most of which were located in the 3'UTR (Supplementary Figure S3D). The RNA editing profiles of most viral pneumonia were similar to those in IBP pneumonia, with Adar up-regulated and Adarb1 down-regulated (Supplementary Figures S3A, B). The Venn plots showed 12 editing sites in 6 edited genes and 3544 editing sites in 688 edited genes exclusively present in uninfected and infected lung tissues, respectively (Figures 4A, B). Likewise, the editing level, the number of editing sites, and genes were increased after viral infection (Figures 4C–E). However, Adar in H3N2 infection was down-regulated (Supplementary Figures S3E, F–H), which was consistent with its overall changes in editing sites and levels (Supplementary Figures S3F–H). Notably, numerous shared sites were found among viral pneumonia datasets (Figure 4F). As shown in Figure 4G, two significantly cis-regulatory DRE sites with up-regulated RNA editing levels in the infected group were found in CTP synthase 1 (Ctps) (Ctps: chr4:120540377, Spearman r = 0.79, P = 2.3 × 10⁻¹⁵) and Terminal Nucleotidyltransferase 5C (Tent5c) (Tent5c: chr3:100468475, Spearman r = 0.65, P = 5.0 × 10⁻⁹). Two cis-regulatory DRE sites with down-regulated RNA editing levels were also observed in Fas Associated Via Death Domain (Fadd) (Fadd: chr7:144579646, Spearman r = - 0.71, P = 9.1× 10⁻⁸) and Lysine Methyltransferase 2D (Kmt2d: chr15:98852368, coefficient r = - 0.47, P = 4× 10⁻³) (Figure 4H). Enrichment analysis and GSEA showed that these DRE genes were mainly involved in the regulation of the triglyceride biosynthetic process, TNF signaling pathway, influenza A, and response to endogenous stimulus (Figures 4I–K).

By comparing IBP and viral pneumonia, we found that RNA editing in viral infections led to a higher proportion of 3′-UTR variants (Figure 5A). And most of the DRE sites were unique to IBP or viral pneumonia (Figures 5B, C). Interestingly, Spearman correlation analysis revealed more DRE sites were correlated with Adar compared to Adarb1 in both IBP and viral infections (Supplementary Figure S4). Some DRE sites changed consistently in terms of editing level between the two types of infections (Figure 5D), which might further regulate gene expression. In line with this, similar up- or down-regulation in gene expression found in these DRE genes suggested mechanisms common to both types of infections (Table 2). Most of the shared DRE genes were hyper-edited. For example, Schlafen 5 (Slfn5) contained six A-to-I RNA editing sites (Slfn5: chr11:82962566, 82963283, 82963634, 82963686, 82962655 and 82962584). The shared GO and KEGG pathways enriched by DRE genes between IBP and viral infections were mainly related to myelocyte-mediated immunity, autophagy, apoptosis, lysosomes, and small GTPases. (Figures 5E-G and Supplementary Tables S6-S9). Therefore, such findings showed RNA editing changes shared by the two types of infections.

RNA editing could be used in the diagnosis of diseases (41, 42). Therefore, we focused on the change of RNA editing profiles that have certain consensus and specificity, which were used for the diagnosis of related infectious disease models. To determine DRE sites with diagnostic significance, we first performed random forest analysis of the identified DRE sites and obtained the top 30 significant DRE sites (Figures 6A, D). We selected those sites that were only present in either IBP or viral pneumonia for receiver operating characteristic curve (ROC) analysis, kept those with area under curve (AUC) > 0.85 for linear regression analysis (Figures 6B, E) and obtained two diagnostic curves for the two types of infections, respectively (Figures 6C, F). Further comparison of these DRE sites between IBP and viral infections obtained seven sites with AUC > 0.85, which were included in a diagnostic model to predict the type of infections (Figures 6G– I). And the results showed that combined analysis of these sites had high sensitivity and specificity in distinguishing IBP infections from viral infections.

Raw data of RNA-Seq were downloaded from the European Nucleotide Archive (ENA) of the European Molecular Biology Laboratory (https://www.ebi.ac.uk/ena). The details of all bacterial and viral infection datasets can be accessed in Table 1.

The process of RNA sequencing reads was conducted as previously described (73). In brief, the raw sequencing data analyzed by FastQC for quality control were aligned and mapped to the mouse genome (UCSC mm10) using RNA STAR (version 2.7.0e) (74). SamTools (version 1.16) was used to filter the reads by removing optic duplications (75), and only reads uniquely mapped were kept. Base quality score recalibration was then performed with the resulting BAM files by using GATK (version 4.1.3) and following the best practice workflows recommended by the documentation (76).

Single nucleotide variants (SNV) were called by using VarScan (version 2.4.4) (77). The variant calling criteria were set as follows: base quality ≥25, total sequencing depth ≥10, alternative allele depth ≥2, and alternative allele frequency (AAF) ≥1%, and possible false positive SNVs were filtered and removed using VarScan with default parameters. SNVs were annotated using the Ensembl Variant Effect Predictor (VEP) (78). SNVs were further filtered and removed according to criteria described previously (73): (1) located in homopolymer runs ≥ 5 nucleotides (nt), simple repeats, in the mitochondria, within 6 nt from splice junctions, within 1 nt from insertions or deletions, or within 4% to the ends of reads; (2) annotated in the dbSNP database Build 142 unless annotated as RNA editing sites in the REDIportal V2.0 database (79) (3); more than 90% of all samples had an AAF equal to 100% or between 40% and 60% (80). High-confidence A-to-I (G) RNA editing SNVs (including A-to-G genomic SNVs on the coding strand and T-to-C genomic SNVs on the opposite strand) were defined either as known RNA editing sites in the REDIportal V2.0 database, or located in genic regions and detected in at least 2 samples with editing levels ≥1%.

Pseudo-counts of gene expression were calculated from the RNA-Seq alignment files using FeatureCounts (81), and transcripts read per thousand bases per million mappings (TPM) were then obtained for each gene using edgeR (version 3.7) (82).

Principal component analysis (PCA) was performed using the function prcomp in R (version 4.2.1) and visualized using the ggplot2 package (version 2.2.1). Heatmaps were plotted using the Pheatmap package in R (version 4.2.1).

Random Forest (83) was used to identify RNA editing sites as biomarkers with high sensitivity and specificity for the diagnosis of infection types. The receiver operator characteristic curve analysis was performed and the area under the curve (AUC) was calculated using the pROC package of R (84).

The enrichment analysis of genes with RNA editing including gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathways were conducted using Enrichr (85).

The GSEA version 4.2.3 software and dataset were used to function of genes based on the GSEA website MSIGDB database (https://www.gsea-msigdb.org/gsea/msigdb/mouse_geneset_resources.jsp) (86), using a default weighted enrichment method with 1000 permutations. Enrichment with false discovery rate (FDR) < 0.25, nominal P-value < 0.05, and |normalized enrichment score (NES)|> 1 were considered significant. NES indicated the analysis results across gene sets. Pairwise P-values were calculated using the non-parametric Kruskal-Wilcoxon test followed by the Tukey post-hoc test.

The generalized linear model (GLM) method and likelihood ratio test were used to compare the intergroup RNA editing levels. The Student’s t-test was used to compare gene expression levels. The Spearman correlation was used to analyze the correlation between RNA editing and gene expression, and correlation coefficients (r) and P-values were calculated. The statistical significance level was set at P < 0.05.