key: cord-0807086-80zq4j7j authors: Di Giorgio, Salvatore; Martignano, Filippo; Torcia, Maria Gabriella; Mattiuz, Giorgio; Conticello, Silvestro G. title: Evidence for host-dependent RNA editing in the transcriptome of SARS-CoV-2 date: 2020-06-17 journal: Sci Adv DOI: 10.1126/sciadv.abb5813 sha: 0f3e6c33db21c0b9851819d053c8cc81d585515b doc_id: 807086 cord_uid: 80zq4j7j 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 [singlestranded 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) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) . On the other hand, APOBECs target the viral genome, typically DNA intermediates (19) (20) (21) (22) (23) (24) (25) (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) (12) (13) (14) (15) (30) (31) (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 singlenucleotide variants (SNVs) on these eight samples (33, 34) using REDItools 2 (35) (36) (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) (45) (46) (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 (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) (53) (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) (20) (21) (22) (23) (24) (25) (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 nonsensemediated 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 import ant 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. 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: 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. 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: 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: 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. 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) . 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. Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/6/25/eabb5813/DC1 A new coronavirus associated with human respiratory disease in China The epitranscriptome and innate immunity A-to-I RNA editing-Immune protector and transcriptome diversifier APOBECs and virus restriction. Virology 479-480 Modeling the embrace of a mutator: APOBEC selection of nucleic acid ligands New antiviral pathway that mediates hepatitis C virus replicon interferon sensitivity through ADAR1 A-to-G hypermutation in the genome of lymphocytic choriomeningitis virus Double-stranded RNA adenosine deaminases enhance expression of human immunodeficiency virus type 1 proteins Evidence for ADAR-induced hypermutation of the Drosophila sigma virus Editing of HIV-1 RNA by the double-stranded RNA deaminase ADAR1 stimulates viral infection Tipping the balance: Antagonism of PKR kinase and ADAR1 deaminase functions by virus gene products Multiple levels of PKR inhibition during HIV-1 replication ADAR2 editing enzyme is a novel human immunodeficiency virus-1 proviral factor Protein kinase PKR and RNA adenosine deaminase ADAR1: New roles for old players as modulators of the interferon response ADARs: Viruses and innate immunity ADARs and the balance game between virus infection and innate immune cell response A-to-I editing of Malacoherpesviridae RNAs supports the antiviral role of ADAR1 in mollusks Selection, recombination, and G----A hypermutation of human immunodeficiency virus type 1 genomes DNA deamination mediates innate immunity to retroviral infection Extensive editing of a small fraction of human T-cell leukemia virus type 1 genomes by four APOBEC3 cytidine deaminases G to A hypermutation of hepatitis B virus Evidence for editing of human papillomavirus DNA by APOBEC3 in benign and precancerous lesions Genetic editing of herpes simplex virus 1 and Epstein-Barr herpesvirus genomes by human APOBEC3 cytidine deaminases in culture and in vivo Extremely high mutation rate of HIV-1 in vivo Characterization of BK polyomaviruses from kidney transplant recipients suggests a role for APOBEC3 in driving in-host virus evolution Antiviral function of APOBEC3G can be dissociated from cytidine deaminase activity Deep sequencing of HIV-1 reverse transcripts reveals the multifaceted antiviral functions of APOBEC3G APOBEC3-mediated restriction of RNA virus replication Long-term restriction by APOBEC3F selects human immunodeficiency virus type 1 variants with restored Vif function APOBEC3G contributes to HIV-1 variation through sublethal mutagenesis Infectious vaccine-derived rubella viruses emerge, persist, and evolve in cutaneous granulomas of children with primary immunodeficiencies RNA based mNGS approach identifies a novel human coronavirus from two individual pneumonia cases in 2019 Wuhan outbreak Genomic diversity of SARS-CoV-2 in coronavirus disease 2019 patients REDItools: High-throughput RNA editing detection made easy Investigating RNA editing in deep transcriptome datasets with REDItools and REDIportal HPC-REDItools: A novel HPC-aware tool for improved large scale RNA-editing analysis JACUSA: Site-specific identification of RNA editing events from replicate sequencing data Infidelity of SARS-CoV Nsp14-exonuclease mutant virus replication is revealed by complete genome sequencing Coronaviruses lacking exoribonuclease activity are susceptible to lethal mutagenesis: Evidence for proofreading and potential therapeutics A-to-I RNA editing occurs at over a hundred million genomic sites, located in a majority of human genes RNA-DNA differences in human mitochondria restore ancestral form of 16S ribosomal RNA Double-stranded RNA adenosine deaminases ADAR1 and ADAR2 have overlapping specificities Genome-wide identification of human RNA editing sites by parallel DNA capturing and sequencing Accurate identification of A-to-I RNA editing in human by transcriptome sequencing Profiling RNA editing in human tissues: Towards the inosinome Atlas REDIportal: A comprehensive database of A-to-I RNA editing events in humans Transcriptome-wide sequencing reveals numerous APOBEC1 mRNA-editing targets in transcript 3' UTRs RNA editors, cofactors, and mRNA targets: An overview of the C-to-U RNA editing machinery and its implication in human disease Genome-wide quantification of ADAR adenosineto-inosine RNA editing activity Widespread cleavage of A-to-I hyperediting substrates Hyperediting of human T-cell leukemia virus type 2 and simian T-cell leukemia virus type 3 by the dsRNA adenosine deaminase ADAR-1 A genome-wide map of hyper-edited RNA reveals numerous new sites The apolipoprotein B mRNA editing complex performs a multifunctional cycle and suppresses nonsense-mediated decay Transient overexpression of exogenous APOBEC3A causes C-to-U RNA editing of thousands of genes Population stratification of a common APOBEC gene deletion polymorphism A common deletion in the APOBEC3 genes and breast cancer risk Trimmomatic: A flexible trimmer for Illumina sequence data Fast and accurate short read alignment with Burrows-Wheeler transform Genome Project Data Processing Subgroup, The sequence alignment/ map format and SAMtools Advanced multi-sample quality control for high-throughput sequencing data Software for computing and annotating genomic ranges MUSCLE: Multiple sequence alignment with high accuracy and high throughput 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. 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.