The COVID-19 outbreak has become a global health risk, and understanding the response of the host to the SARS-CoV-2 virus will help to combat the disease. RNA editing by host deaminases is an innate restriction process to counter virus infection, but it is not yet known whether this process operates against coronaviruses. Here, we analyze RNA sequences from bronchoalveolar lavage fluids obtained from coronavirus-infected patients. We identify nucleotide changes that may be signatures of RNA editing: adenosine-to-inosine changes from ADAR deaminases and cytosine-to-uracil changes from APOBEC deaminases. Mutational analysis of genomes from different strains of Coronaviridae from human hosts reveals mutational patterns consistent with those observed in the transcriptomic data. However, the reduced ADAR signature in these data raises the possibility that ADARs might be more effective than APOBECs in restricting viral propagation. Our results thus suggest that both APOBECs and ADARs are involved in coronavirus genome editing, a process that may shape the fate of both virus and patient.
Emerging viral infections represent a threat to global health, and the recent outbreak of novel coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2, novel coronavirus, 2019-nCoV) exemplifies the risks (1, 2). As viruses are obligate intracellular parasites, organisms have evolved innate immune mechanisms to sense and counter the viruses. Among these mechanisms, RNA and DNA editing mediated by endogenous deaminases can provide a potent defense against specific viruses. Two deaminase families are present in mammalian species: The ADARs (adenosine deaminases that act on RNA) target double-stranded RNA (dsRNA) for deamination of adenines into inosines (A-to-I) (3, 4), and the APOBECs deaminate cytosines into uracils (C-to-U) on single-stranded nucleic acids [single-stranded DNA (ssDNA) and single-stranded RNA (ssRNA)] (5, 6). During viral infections, ADARs act either directly, through hypermutation of the viral RNA, or indirectly, through editing of host transcripts that modulate the cellular response (7–18). On the other hand, APOBECs target the viral genome, typically DNA intermediates (19–26), either through C-to-U hypermutation or through a non-enzymatic path that interferes with reverse transcription (27, 28). Although some APOBEC3 proteins can interfere in vitro with Coronaviridae replication, it is not clear whether their enzymatic activity is involved (29). Ultimately, though, these restriction systems can also be exploited by the viruses to support their infectivity and increase their evolutionary potential (9, 11–15, 30–32).
To assess whether RNA editing could be involved in human host responses to SARS-CoV-2 infections, we started from publicly available RNA sequencing datasets from bronchoalveolar lavage fluids (BALF) obtained from patients diagnosed with COVID-19. While transcriptomic data for all samples could be aligned to the SARS-CoV-2 reference genome, the quality of the sequencing varied and only eight samples had coverage and error rates suitable for the identification of potentially edited sites (data S1). We called single-nucleotide variants (SNVs) on these eight samples (33, 34) using REDItools 2 (35–37) and JACUSA (38) using the following thresholds: reads supporting the SNV ≥4, allelic fraction ≥0.5%, coverage ≥20, quality of the reads >25, base quality >35 (fig. S1A). The two pipelines gave comparable results with ~50% of the SNV positions called by both (figs. S1B and S2). We identified 910 SNVs common to REDItools 2 and JACUSA, ranging from 24 to 238 SNVs per sample (Fig. 1 and data S3). Given the thresholds used to call the SNV, samples with lower sequencing depths displayed lower numbers of SNVs.
While the weight of each SNV type varies across samples (Fig. 1), a bias toward transitions is always present, which is even more evident when all mutational data are pooled (Fig. 2, A and B). This pattern holds true even when only SNVs recurring in more samples are considered (Fig. 2C).
The SNV frequency and number of transversions are compatible with the mutation rates observed in coronaviruses [10–6/−7; (39)] and commonly associated to the RNA-dependent RNA polymerases (RdRps). RdRps are error prone and are considered the main source of mutations in RNA viruses. However, the coronavirus nsp14-ExoN gene provides a form of error correction (40), which is probably the reason mutation rates in coronaviruses are lower than those observed in RNA viruses with smaller genomes. The mutational spectrum in SARS quasispecies presents a very weak bias toward U-to-G. Inactivation of nsp14-ExoN error correction reveals the mutational spectrum of the RdRp, which is quite different from the pattern we observe (i.e., main changes are C-to-A, followed by U-to-C, G-to-U, A-to-C, and U-to-G) (41). Hence, we would consider that SNVs deriving from RdRp errors represent a marginal fraction of the SNVs in the SARS-CoV-2 samples.
The bias toward transitions—mainly A>G/T>C changes—resembles the pattern of SNVs observed in human transcriptomes (42) or in viruses (8, 10, 18), where A>G changes derive from deamination of A-to-I mediated by the ADARs. It is thus likely that the A>G/T>C changes seen in SARS-CoV-2 are also due to the action of ADARs.
C>T and G>A SNVs are the second main group of changes and could derive from APOBEC-mediated C-to-U deamination. Unlike A-to-I editing, C-to-U editing is a relatively rare phenomenon in the human transcriptome (42), and with regard to viruses, it has been associated only with positive-sense ssRNA rubella virus (32), where C>T changes represent the predominant SNV type. The observation that only A-to-I editing is present in RNA viruses that infect nonvertebrate animals, where RNA-targeting APOBECs are not present (10, 18), supports the hypothesis that APOBECs are involved in the RNA editing of this human-targeting virus.
A third group of SNVs, A>T/T>A transversions, is also present in these samples. While this type of SNV has been reported in other genomic studies (43), its origin is still unknown.
A>G and T>C changes are evenly represented with respect to SNV frequency (Fig. 2A), the number of unique SNVs (Fig. 2, B and C), and their distribution across the viral genome (Fig. 2D). As ADARs target dsRNA, this suggests that dsRNA encompasses the entire genome. While dsRNA in human transcripts is often driven by inverted repeats, the most likely source of dsRNA in the viral transcripts is replication, where both positive and negative strands are present and can result in wide regions of dsRNA.
Unlike A-to-I changes, C-to-U changes are biased toward the positive-sense strand (Fig. 2, B to D; P < 0.0001). Because ADARs and APOBECs selectively target dsRNA and ssRNA, this distribution could arise from the presence at all times of RNA in a dynamic equilibrium between double-strandedness—when negative-sense RNA is being transcribed—and single-strandedness—when nascent RNA is released. Although some areas seem to bear fewer SNVs, these reduced SNV frequencies might be related to lower sequencing depth in those regions.
As APOBEC deaminases preferentially target cytosines within specific sequence contexts, we analyzed the nucleotide context of A-to-I and C-to-U SNVs in the viral genome (Fig. 3, A and B). A slight depletion of G bases in position −1 is present at A-to-I edited positions. This depletion is not as strong as the signal previously reported in human transcripts (44–47). The low editing frequencies we observe resembles the editing present on human transcripts containing Alu sequences, which were found in a limited number in those early datasets. There is no evidence of a sequence context preference if we use a larger dataset such as REDIportal (48), which includes >1.5 M sites in Alu repeats (fig. S3).
On the other hand, C-to-U changes preferentially occur downstream from uridines and adenosines, within a sequence context that resembles the one observed for APOBEC1-mediated deamination ([AU]C[AU]) (49, 50).
We then aligned available genomes from SARS-CoV-2, Middle-East respiratory syndrome–related coronavirus (MERS-CoV), and SARS-CoV to test whether RNA editing could be responsible for some of the mutations acquired through evolution. The genomic alignments reveal that a substantial fraction of the mutations in all strains could derive from enzymatic deaminations (Fig. 4, A to C), with a prevalence of C-to-U mutations, and a sequence context compatible with APOBEC-mediated editing also exists in the genomic C-to-U SNVs (Fig. 4, D to F).
Our data source—metagenomic sequencing—raises the question whether the low-level editing we observe (~1%) reflects the actual levels of editing of viral transcripts within human cells. Aside from a small fraction of cellular transcripts edited at high frequency, most ADAR-edited sites in the human transcriptome (typically inside Alu sequences) present editing levels of ~1% (4, 42, 51). It has been shown that a fraction of the cellular transcripts are hyperedited by ADARs (52–54). While we were unable to observe hyperedited reads in the metagenomic samples, it is possible that hyperedited transcripts fail to be packaged into the virus.
With regard to APOBEC-mediated RNA editing, its detection in the viral transcriptomes is already indicative, as this type of editing is almost undetectable in human tissues (42). This enrichment points either toward an induction of APOBECs triggered by coronavirus infection or to specific targeting of the APOBECs to the viral transcripts. APOBECs have been proved effective against many viral species in experimental conditions, yet, until now, their mutational activity in clinical settings has been shown only in a handful of viral infections (19–26) through DNA editing and, in rubella virus, on RNA (32).
As in rubella virus, we observe a bias in APOBEC editing toward the positive-sense strand. This bias and the low editing frequencies might be indicative of the dynamics of the virus, from transcription to selection of viable genomes. It is reasonable to assume that sites edited on the negative-sense strand will result in a mid-level editing frequency, as not all negative-sense transcripts will be edited (Fig. 5A). On the other hand, editing of the positive-sense strand can occur upon entry of the viral genome, thus yielding high-frequency editing (Fig. 5B), or after viral genome replication, resulting in low-frequency editing (Fig. 5C). The lack of a sizable fraction of highly edited C>T SNVs suggests that APOBEC editing occurs late in the viral life cycle (Fig. 5C). Yet, because they occur earlier, G>A SNVs should be closer in number to C>T SNVs and with higher levels of editing, which is not what we observe (Fig. 2, A to C). The overrepresentation of C>T SNVs could be due to an imbalance toward positive-sense transcripts, as these are continuously generated from the negative-sense ones (and double-stranded hybrid RNAs are lost). However, the editing frequencies of G>A SNVs should be much higher, as G>A SNVs are generated upstream to the C>T ones. A more fitting explanation is that editing of the negative-sense transcripts results in loss of the edited transcript (Fig. 5D), possibly because editing triggers nonsense-mediated decay (55), thus lowering the chances of the edited site to be transmitted.
Because most of the APOBECs are unable to target RNA, the only well-characterized cytidine-targeting deaminases are APOBEC1, mainly expressed in the gastrointestinal tract, and APOBEC3A (56), whose physiological role is not clear. As with A-to-I editing, it will be important to assess the true extent of APOBEC RNA editing in infected cells.
The functional meaning of RNA editing in SARS-CoV-2 is yet to be understood: In other contexts, editing of the viral genome determines its demise or fuels its evolution. For DNA viruses, the selection is indirect, as genomes evolve to reduce potentially harmful editable sites [e.g., (18)], but for RNA viruses, this pressure is even stronger, as RNA editing directly affects the genetic information and efficiently edited sites disappear.
A comparison of the SNV datasets from the transcriptomic and genomic analyses reveals a different weight of A-to-I and C-to-U changes (Figs. 2B and 4A), with an underrepresentation of A-to-I in the viral genomes. As our analysis underestimates the amount of editing due to the strict parameters used, the underrepresentation of A-to-I changes could be explained by the possibility that A-to-I editing is more effective in restricting viral propagation, thus reducing the number of viral progeny showing evidence of these changes. In contrast, the remnants of less effective C-to-U editing are retained in viral progeny and get fixed during viral adaptation.
An analysis of mutation outcomes is difficult due to the low numbers of events collected so far, but there are some possibly suggestive trends (data S2). C-to-U changes leading to stop codons are overrepresented in the transcriptomic data but—as expected—disappear in the genomic dataset. This might point—again—to an antiviral role for these editing enzymes. There is also an underrepresentation of C>T missense mutations, but its meaning is difficult to interpret.
Last, this analysis is a first step in understanding the involvement of RNA editing in viral replication, and it could lead to clinically relevant outcomes: (i) If these enzymes are relevant in the host response to coronavirus infection, a deletion polymorphism quite common in the Chinese population, encompassing the end of APOBEC3A and most of APOBEC3B (57, 58), could play a role in the spread of the infection. (ii) Because RNA editing and selection act orthogonally in the evolution of the viruses, comparing genomic sites that are edited with those that are mutated could lead to the selection of viral regions potentially exploitable for therapeutic uses.
MATERIALS AND METHODS
RNA sequencing data available from projects PRJNA601736, PRJNA603194, and PRJNA605907 were downloaded from the National Center for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/sra/) using the FASTQ-dump utilities from the SRA-toolkit with the following command line:
prefetch -v SRR* && fastq-dump --outdir /path_dir/ | --split-files /path_dir/SRR*.sra
Because most of the reads of samples from PRJNA605907 were missing their mate, forward-reads and reverse-reads from these samples have been merged in a single FASTQ, which is treated as a single-end experiment. Details of the sequencing runs are summarized in data S1.
SRR11059940, SRR11059941, SRR11059942, and SRR11059945 showed a reduced quality of the sequencing in the terminal part of the reads. We used TRIMMOMATIC (59) to trim the reads of those samples to 100 base pairs (bp), with the following command line:
rimmomatic SE SRR*.fastq SRR*.trimmed.fastq CROP:100
We aligned the FASTQ files using Burrows-Wheeler Aligner (60) using the official sequence of SARS-CoV-2 (NC_045512. 2) as reference genome. After the alignments, BAM files were sorted using SAMtools (61).
The command line used for paired-end samples is as follows:
bwa mem NC_045512.2.fa SRR*_1.fastq SRR*_2.fastq | samtools sort –O BAM -o SRR*_.bam
The command line used for single-end samples is as follows:
bwa mem NC_045512.2.fa SRR*.fastq | samtools sort –O BAM -o SRR*_.bam
The aligned bams have been analyzed with QUALIMAP (62). Because of a high error rate reported by QUALIMAP, samples SRR11059943 and SRR10971381 have been removed from the analysis.
python2.7 reditools.py -f SRR*.bam -o SRR10903401_stat_table_allPos.txt -S -s 0 -os 4 -m /homol_site/SRR*_homopol.txt -c SRR*_homopol.txt -r /Reference/NC_045512.2.fa -a SRR*_stat_table_allPos.txt -q 25 -bq 35 -mbp 15 -Mbp 15
jacusa call-1 -p 20 -r SRR*.vcf -a B,I,Y -s -f V -q 35 -m 25 SRR*.srt.bam
With regard to REDItools 2, we removed all SNVs within 15 nucleotides from the beginning or the end of the reads to avoid artifacts due to misalignments.
To avoid potential artifacts due to strand bias, we used the AS_StrandOddsRatio parameter, calculated following GATK guidelines (https://gatk.broadinstitute.org/hc/en-us/articles/360040507111-AS-StrandOddsRatio), and any mutation with an AS_StrandOddsRatio > 4 has been removed from the dataset.
Bcftools (61) has been used to calculate total allelic depths on the forward and reverse strand (ADF and ADR) for AS_StrandOddsRatio calculation, with the following command line:
mpileup -a FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP -O v -A -C -I -d 1000000 -q 25 -Q 35 -f NC_045512.2.fa -o SRR*.vcf SRR*.srt.bam
Mutations common to the datasets generated by REDItools 2 and JACUSA were considered (n = 910; fig. S2 and data S3). The threshold we used to filter the SNVs is based on minimum coverage (20 reads), number of supporting reads (at least four mutated reads), allelic fraction (0.5%), quality of the mapped reads (>25), and base quality (>35). In the dataset, there were only six SNVs with allelic fractions in the range of 30 to 85% (C>T, 1; T>C, 3; G>T, 2). Because there were no SNVs with higher allelic fractions, we presume that all samples originated from the same viral strain.
Recurring SNVs have been defined as the SNVs present in at least two samples. To overcome the problem of samples with lower sequencing depth, we used the positions of the SNVs common to both REDItools 2 and JACUSA to call again the SNVs irrespectively of the number of supporting reads.
R packages (Biostrings, rsamtools, ggseqlogo ggplot2, and splitstackshape) and custom Perl scripts were used to handle the data.
Sequence context analysis
Logo alignments were calculated using ggseqlogo, using either the pooled dataset or the dataset of recurring SNVs. Logo alignments of the human edited sites were performed using ADAR sites from REDIportal (48) that were shared by at least four samples. SARS-CoV-2, SARS, and MERS genomic data were prepared for the Logi alignment using the GenomicRanges R package (63).
SNV calling in genomic data from SARS-CoV-2, SARS, and MERS
The viral genomic sequences of MERS (taxid:1335626) and SARS (taxid:694009) were selected on NCBI Virus (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/) using the following query: Host : Homo Sapiens (human), taxid:9606; -Nucleotide Sequence Type: Complete. They were aligned using the “Align” utility. Consensus sequences of SARS and MERS genomes were built using the “cons” tool from the EMBOSS suite (http://bioinfo.nhri.org.tw/gui/) with default settings. SARS-CoV-2 genomic sequences were downloaded from GISAID (https://www.gisaid.org/) and aligned with MUSCLE (64).
SNVs have been called with a custom R script, by comparing viral genome sequences to the respective consensus sequence or, for SARS-CoV-2, to the NC_045512.2 reference sequence. SNVs, viral consensus sequences, and Coronaviridae genome sequence identifiers are provided in data S3 to S5.
SNVs (from both genomic and somatic SNV sets) occurring on coding sequences have been annotated with custom R scripts to determine the outcome of the nucleotide change (nonsense/missense/synonymous mutation). A summary is reported in data S2.
fisher.test() function from the R base package has been used for all the statistical tests. To test the significance of C-to-U bias on the positive strand, we compared C>T/G>A SNV counts to the count of C/G bases on the reference genome. For P values of “RNA vs Reference,” “DNA vs Reference,” and “genome vs RNA,” 2 × 2 contingency tables have been generated as shown in data S2.
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).