Research ArticleCORONAVIRUS

Evidence for host-dependent RNA editing in the transcriptome of SARS-CoV-2

See allHide authors and affiliations

Science Advances  17 Jun 2020:
Vol. 6, no. 25, eabb5813
DOI: 10.1126/sciadv.abb5813

Abstract

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.

INTRODUCTION

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 (718). On the other hand, APOBECs target the viral genome, typically DNA intermediates (1926), 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, 1115, 3032).

RESULTS

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 (3537) 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.

Fig. 1 SNVs identified in SARS-CoV-2 transcriptomes.

The bar charts show the number of SNVs identified in each SARS-CoV-2 transcriptome for each SNV type (e.g., A>C, AC). The sequencing depth for each sample is indicated.

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).

Fig. 2 SNV identified in SARS-CoV-2 transcriptomes.

(A) Allelic fraction and (B) number of SNVs for each nucleotide change in the entire dataset and (C) for SNVs recurring in at least two samples. (D) Distribution of SNVs across the SARS-CoV-2 genome. A-to-G (blue) and C-to-U (red) SNVs are grouped in 400-nucleotide (nt) bins and plotted above (AG and CT) or below the line (TC and GA) based on the edited strand. Dots (white/black) indicate recurring SNVs. Genetic organization of SARS-CoV-2 (top). The dark/white shading indicates the viral coding sequences; coverage distribution of all analyzed samples (bottom).

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 (4447). 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).

Fig. 3 Sequence contexts for SARS-CoV-2 RNA edited sites.

(A) Local sequence context for A-to-I and C-to-U edited sites in the viral transcriptome and (B) for recurring sites.

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).

Fig. 4 Nucleotide changes across Coronaviridae strains.

(A to C) Number of SNVs for each nucleotide change and (D to F) local sequence context for C-to-U edited sites in genome alignments from SARS-CoV-2 (A and D), human-hosted MERS-CoV (B and E), and human-hosted SARS-CoV (C and F).

DISCUSSION

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 (5254). 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 (1926) 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.

Fig. 5 Model of APOBEC RNA editing on SARS-CoV-2 transcriptome.

The four panels model the editing frequencies and the C>U/G/A ratios expected from four different scenarios: (A) C-to-U editing on the negative-sense transcripts, (B) “early” editing on the viral genomes before viral replication, (C) “late” editing after viral replication, and (D) “late” editing after viral replication with loss of negative-sense transcripts. Red dots indicate editing on the positive-sense transcript; orange dots indicate editing on the positive-sense transcript. Green and blue segments indicate positive- and negative-sense viral transcripts, respectively.

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

Sequencing data

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.

Data preprocessing

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.

SNV calling

A diagram of the entire pipeline is shown in fig. S1A. We used REDItools 2 (35, 37) and JACUSA (38) to call the SNVs using the following command line:

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.

Data manipulation

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.

SNV annotation

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.

Statistical analysis

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.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/25/eabb5813/DC1

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.

REFERENCES AND NOTES

Acknowledgments: In memory of Li Wenliang, Carlo Urbani, and all the doctors and health workers who endangered their lives in the fight against epidemics. We acknowledge and thank all authors that are sharing their data. We gratefully acknowledge the authors, originating and submitting laboratories, of the sequences from GISAID’s EpiFlu Database on which this research is based. The list is detailed in data S6. We thank Ernesto Picardi for the support with REDItools2. Funding: The research was supported by grants from Ministero della Salute (PE-2013-02357669) and from AIRC (IG-17701). Author contributions: Conceptualization: S.D.G., F.M., M.G.T., G.M., and S.G.C.; formal analysis, investigation, and software: S.D.G. and F.M.; visualization: G.M.; writing (original draft): G.M. and S.G.C.; writing (review and editing): S.D.G., F.M., M.G.T., G.M., and S.G.C. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Sequencing and genomic data are available through NCBI SRA and Virus repositories and GISAID. The SNVs dataset is also available through the UCSC genome browser at http://genome.ucsc.edu/s/Max/conticello2020.
View Abstract

Stay Connected to Science Advances

Navigate This Article