Emerging Microbes & Infections 2021, VOL. 10 https://doi.org/10.1080/22221751.2021.1969868 Host-dependent editing of SARS-CoV-2 in COVID-19 patients Josep Gregori a,b,c*, Maria Francesca Cortese d*, Maria Piñana e, Carolina Campos a,b, Damir Garcia-Cehic a,b, Cristina Andrés e, Josep Francesc Abril f, Maria Gema Codina g, Ariadna Rando g, Juliana Esperalba g, Elena Sulleiro g, Joan Joseph g, Narcís Saubí d, Sergi Colomer-Castell a, Mari Carmen Martin e, Carla Castilloe, Juan Ignacio Esteban a, Tomas Pumarola g,h, Francisco Rodriguez-Frias b,d,h, Andrés Antón e,h and Josep Quer a,b a Liver Diseases-Viral Hepatitis, Liver Unit, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Barcelona Hospital Campus, Barcelona, Spain; bCentro de Investigación Biomédica en Red de Enfermedades Hepáticas y Digestivas (CIBERehd), Instituto de Salud Carlos III, Madrid, Spain; cRoche Diagnostics SL, Barcelona, Spain; dBiochemistry and Microbiology Departments, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Barcelona Hospital Campus, Barcelona, Spain; eRespiratory Viruses Unit, Microbiology Department, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Barcelona Hospital Campus, Barcelona, Spain; fComputational Genomics Lab, Genetics, Microbiology and Statistics Department, Institut de Biomedicina (IBUB), Universitat de Barcelona, Barcelona, Spain; gMicrobiology Department, Vall d’Hebron Hospital Universitari, Vall d’Hebron Barcelona Hospital Campus, Barcelona, Spain; hUniversitat Autònoma de Barcelona, Bellaterra, Spain ABSTRACT A common trait among RNA viruses is their high capability to acquire genetic variability due to viral and host mechanisms. Next-generation sequencing (NGS) analysis enables the deep study of the viral quasispecies in samples from infected individuals. In this study, the viral quasispecies complexity and single nucleotide polymorphisms of the SARS-CoV-2 spike gene of coronavirus disease 2019 (COVID-19) patients with mild or severe disease were investigated using next-generation sequencing (Illumina platform). SARS-CoV-2 spike variability was higher in patients with long- lasting infection. Most substitutions found were present at frequencies lower than 1%, and had an A → G or T → C pattern, consistent with variants caused by adenosine deaminase acting on RNA-1 (ADAR1). ADAR1 affected a small fraction of replicating genomes, but produced multiple, mainly non-synonymous mutations. ADAR1 editing during replication rather than the RNA-dependent RNA polymerase (nsp12) was the predominant mechanism generating SARS-CoV-2 genetic variability. However, the mutations produced are not fixed in the infected human population, suggesting that ADAR1 may have an antiviral role, whereas nsp12-induced mutations occurring in patients with high viremia and persistent infection are the main source of new SARS-CoV-2 variants. ARTICLE HISTORY Received 6 April 2021; Revised 30 July 2021; Accepted 12 August 2021 KEYWORDS ADAR1; SARS-CoV-2; editing; mutations; quasispecies Introduction recombination events) have been described in several coronaviruses, including SARS-CoV-2 [4–10]. Severe acute respiratory syndrome coronavirus 2 Animal cells have a first-line innate immunity that (SARS-CoV-2), the causal agent of the acute respiratory enables a prompt response to danger signals and syndrome known as coronavirus disease 2019 (COVID- pathogens. Interactions between the specialized cellu- 19), replicates by using its own RNA-dependent RNA lar receptors known as pattern recognition receptors polymerase (nsp12). SARS-CoV-2 expresses an acces- (PRRs) and conserved pathogen-associated molecular sory non-structural protein, nsp14, which has 3’-5’ exo- patterns (PAMPs) rapidly trigger an innate intracellu- nuclease proofreading activity. This feature, which lar antiviral response, usually based on immediate notably limits the acquisition of mutations during repli- activation of the interferon and NF-kb pathways cation of the viral genome, is seen in other corona- [11]. Unlike bacteria, viral DNA and RNA have cer- viruses, such as SARS-CoV and mouse hepatitis virus tain features (unusual in the cellular host) that can [1, 2]. Phylogenetic analysis of SARS-CoV-2 consensus directly activate specific PRRs. One such PRR is ade- sequences has shown only small differences [3], but nosine deaminase acting on RNA1 (ADAR1). other sources of variation (insertions, deletions, and ADAR1 directly attacks the virus by deaminating CONTACT Josep Quer josep.quer@vhir.org, josep.quer@ciberehd.org Virus Translational Research Unit, Liver Diseases, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Barcelona Hospital Campus, Passeig Vall d’Hebron 119-129, Barcelona 08035, Spain; Andrés Antón aanton@vhebron.net Respiratory Viruses Unit, Microbiology Department, Vall d’Hebron Institut de Recerca (VHIR), Vall d’Hebron Hospital Universitari, Pg Vall d’Hebron 119- 129, Barcelona 08035, Spain *These two authors have contributed equally to the manuscript. Supplemental data for this article can be accessed at https://doi.org/10.1080/22221751.2021.1969868. © 2021 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrest- ricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1778 J. Gregori et al. adenosine (A) in double-stranded RNA, thereby pro- contribution of the intracellular mutagenic enzyme, ducing inosine (I) which is misinterpreted as guano- ADAR1, to the variability observed. sine (G) by the splicing and translation machinery. ADAR1 editing of double-stranded RNA (dsRNA) can be specific (deamination of selected A residues) Patients and methods or almost random. In long dsRNA (>100 bp), this Patients editing can affect 50% of all adenosine residues and cause a hypermutation, as has been reported in some The study included samples from 18 patients with single-stranded RNA viruses during replication and mild or severe COVID-19 admitted to the emergency persistent infection [12]. In contrast, short dsRNSs room of Hospital Universitari Vall d’Hebron (∼20–30 bp) and long dsRNAs with bulges, helixes, (HUVH), Barcelona (Spain) in March 2020. SARS- hairpins, and loops are edited selectively and only a CoV-2 infection was diagnosed in the HUVH Micro- few adenosines are specifically chosen [13]. biology Department using two tests, an in-house PCR RNA editing by ADAR1 has been reported in sev- assay using the primer/probe set from the CDC 2019- eral viruses [14], where it may contribute to inhibiting nCoV real-time RT-PCR diagnostic panel (Qiagen, viral infection [15]. However, viruses can also take Germany) and the commercial Allplex 2019-nCoV advantage of the variability introduced by cellular real-time RT-PCR assay (Seegene, Korea). Samples mutagenic enzymes, including ADAR1, to support had been collected from the upper respiratory tract viral infectivity and evolution [16–18]. Adenosine dea- (nasopharyngeal/oropharyngeal swabs or nasophar- mination in codon 196 of the hepatitis delta virus yngeal aspirates). Ten patients had mild symptoms (HDV) antigenome enables the expression of large (absence of influenza-like illness or hypoxia and no delta antigen, which has a key role in viral replication hospitalization requirement), and 8 patients had [19]. In human immunodeficiency virus 1 (HIV-1), severe disease (e.g. pneumonia, hypoxemic respiratory measles virus (MV), vesicular stomatitis virus (VSV), failure, sepsis, cardiomyopathy and arrhythmia, acute and hepatitis D virus (HDV), ADAR1 enhances repli- kidney dysfunction) and were admitted to the inten- cation during the acute phase of infection [20]. It acts sive care unit. All patients had a favourable outcome, by editing the viral substrates – adding A to G substi- with resolution of the infection (Table 1). No viral tutions to the RNA genome [21] – or by inhibiting or bacterial co-infections were reported. dsRNA-dependent protein kinase (PKR) [20]. Due The patients included had no comorbidities other to the large genome size of SARS-CoV-2 (≈30 kb) than acute respiratory syndrome and were all SARS- and its complex tridimensional structure [22], the CoV-2 monoinfected. All cases included occurred ADAR enzyme together with other deaminases may during the first wave of the COVID-19 pandemic, in be involved in the genome editing process [23]. which treatment of severely ill patients consisted SARS-CoV-2 has shown low variability to date, but mainly of respiratory maintenance because there had deep analysis of the viral genome might detect variants been no previous experience in the management of overlooked by traditional sequencing. We hypothesized these patients. that in-depth analysis of the SARS-CoV-2 spike gene in samples obtained from COVID-19 patients at the begin- Methods ning of the pandemic could reveal the underlying mech- anisms causing variants in this part of the genome before SARS-CoV-2 sequencing and quality control they become fixed in the viral quasispecies. Are the sub- Viruses were inactivated by mixing respiratory speci- stitutions, present at low levels, due to replication errors mens with AVL buffer (Qiagen, Hilden, Germany). produced by the viral polymerase or to interactions with Viral RNA was extracted using the QIAmp Viral enzymes inherent to the patients’ immune system? In RNA Mini Kit (Qiagen, Hilden, Germany) following addition, examination of the viral quasispecies in patients with different clinical profiles (mild or severe Table 1. Clinical and demographic characteristics of COVID-19 disease) could reveal clinical parallels. Finally, after deter- patients included in the study. mining the mechanisms that generate replication errors, MILD (n = 10) Severe (n = 8) it could be of interest to compare the mutations found Sample type with those recurrently detected worldwide after several Nasopharyngeal aspirate 30% (3/10) 12.5% (1/8) Nasopharyngeal/oropharyngeal 70% (7/10) 87.5% (7/8) months of the pandemic to study the impact of these SWAB mechanisms on the profile of fixed mutations. RT-PCR threshold, Ct median[IQR]* 23.4 [19.6; 30.8 [24.2; 25.5]* 32.7]* In this study, we analyzed the SARS-CoV-2 spike Sex, % men 30% (3/10) 50% (4/8) gene by next-generation sequencing (NGS) in samples Age, years 42 [28; 50]* 48 [45; 50]* Days in intensive care unit Not admitted 16 [7; 25]* obtained from COVID-19 patients. Our aim was to The asterisk (*) indicates the minimum and maximum value of the inter- meticulously determine the presence of nucleotide vals included in each square bracket. Ct, cycle threshold; IQR, interquar- mutations in the viral genome and investigate the tile range. Emerging Microbes & Infections 1779 the manufacturer’s instructions. The complete spike and posterior amplicons, thanks to the amplicon over- gene was amplified and sequenced into 13 overlap- lapping. No significant conflicts were observed. ping amplicons, as previously reported by Andres Haplotypes in each sample/amplicon were right et al. [10]. trimmed to the next amplicon and the spike gene The RT-PCR amplification and sequencing method was analyzed haplotype-by-haplotype and site-by- has been described [10]. The final library was loaded in site. A variant was defined as a nucleotide that differed a MiSeq Reagent Kit 600V3 cartridge (Illumina, San from that of the corresponding dominant haplotype Diego, CA) and sequenced using the MiSeq platform (master sequence). (Illumina, San Diego, CA), which, in our experience, To study eventual substitution biases, the frequency has the highest resolution power [24, 25]. To obtain of a variant per site was computed as the number of high-quality haplotypes, the sequencing reads under- haplotypes bearing the variant in the overall study went rigorous quality filters [26]. Reads had to have population, without taking into consideration the hap- a minimum of 20 overlapping bases and a maximum lotype frequency in the viral quasispecies of each of 10% mismatches. Reads that did not fulfil these specific sample/amplicon. The Fisher test was used requirements and those with more than 1% of bases to contrast biases in the base substitutions from the below a Phred score of Q33 (0.9995 accuracy) were expected values, and to determine homogeneity vs excluded. independence in the distribution of pairs of Reads were demultiplexed by matching primers, substitutions. allowing a maximum of three mismatches, and the The above-described analyses were conducted for primer sequences at both ends were trimmed off. both the overall study population and for each specific Identical reads were collapsed to haplotypes, and patient. The frequency of each nucleotide, the number their corresponding frequencies were expressed as of polymorphic sites, and the substitution type (tran- read counts. A fasta file was generated with each sition or transversion) were determined for each pool/primer/strand combination. Reverse haplotypes patient. were reverse-complemented. Haplotypes with more than three gaps were removed. Remaining gaps were Characterization of single nucleotide variants repaired based on the most abundant haplotype (mas- (SNVs) by type and frequency ter) in the corresponding quasispecies. Haplotypes The 3’ end of each amplicon was trimmed to avoid were selected only if they had both a forward and redundancies in the overlaps, and the alignment was reverse strand and were present at >0.1% in each used to count the substitutions that occurred. Regard- strand, providing a robust quality control. less of the abundance of each haplotype in the quasis- To compensate for the differing coverage of the pecies, the number of SNVs was computed with samples, we established a reference coverage of respect to the quasispecies dominant haplotype in 70,000 reads and an abundance threshold of 0.15%. each sample/amplicon. Haplotypes present at ≤0.15% were excluded [27]. SNVs were characterized by substitution type and All computations were done in the R language and by frequency. Variant abundance was expressed as platform [28], and in-house scripts were developed the percentage of each SNV within the particular qua- using Biostrings [29], Ape [30], ShortRead [31], and sispecies (per patient and amplicon). To avoid redun- ggseqlogo [32]. dancy, SNV percentages in overlapping regions of the amplicon were averaged. In addition, SNVs were Quasispecies diversity characterized as the frequency per abundance bin; Quasispecies diversity indices were calculated based that is, the number of SNVs in all the samples (per on the characteristics and frequency of all haplotypes patient and amplicon) with abundances in the follow- observed: number of haplotypes per amplicon, num- ing bins: (0%, 0.2%], (0.2%, 0.5%], (0.5%, 1%], (1%, ber of polymorphic sites per amplicon, frequency of 2%], (2%, 5%], (5%,10%], and (10%,50%]. The lower the dominant haplotype per amplicon, rare haplotype bins would be representative of the underlying mech- load (RHL) at 0.1% and at 1% per amplicon, Shannon anism causing the substitutions, whereas the higher entropy, and mutation frequency [27, 33–36]. bins would incorporate the effects of haplotype fitness and selection. Analysis of variants Minority variants were identified per amplicon, haplo- Location of substitutions in the spike gene type, and site, and were categorized by type and fre- secondary structure quency. Overlapping regions were used as controls To study how SNVs could affect the secondary struc- to compare variants identified in consecutive ampli- ture of SARS-CoV-2 RNA in the spike region, we used cons. Potential bias effects due to the primers were the Wuhan-Hu-1 (MN908927.2) spike gene sequence, controlled by inspecting how the forward and reverse with the D614G (A23403G) substitution in the S primers of each amplicon were seen in the previous protein as the backbone, into which each naturally 1780 J. Gregori et al. Table 2. Statistics of diversity indices regarding the haplotype characteristics and abundance within the corresponding quasispecies. Index Condition Min 1st Q Median Mean 3rd Q Max Depth Severe 89,820 132,475 192,410 214,754 267,141 520,239 Mild 71,306 105,405 127,254 146,090 165,103 492,886 # Hpl Severe 1 2 4 7.82 10.25 32 Mild 1 1 4 5.56 7 24 # PolySites Severe 0 1 3 6.73 9.25 31 Mild 0 0 3 4.58 6.75 23 % Master Severe 73.4 96.1 98.9 96.7 99.6 100 Mild 91.9 97.9 99.3 98.6 100 100 RHL@1% Severe 8.19% 11.78% 13.55% 14.76% 16.04% 28.24% Mild 8.0% 11.0% 12.91% 13.63% 15.27% 22.85% RHL@0.1% Severe 7.64% 10.47% 11.66% 11.72% 12.85% 17.17% Mild 7.23% 10.03% 11.18% 11.34% 12.41% 16.40% Shannon Severe 0 0.0263 0.0740 0.1830 0.2592 11,727 Mild 0 0 0.0512 0.0938 0.1402 0.4859 Mf × 10−6 Severe 0 9.28 30.44 96.65 112.93 659.67 Mild 0 0 24.17 41.79 59.99 264.10 Depth, number of reads per amplicon; # Hpl, number of haplotypes per amplicon; Mf, mutation frequency per amplicon; # % Master, percentage of the most abundant haplotype in the amplicon (master sequence); PolySites, number of polymorphic sites per amplicon; RHL@1%, rare haplotype load at 1%; RHL@0.1%, rare haplotype load at 0.1%; Shannon, Shannon entropy. observed nucleotide substitution was inserted. A mini- detected in mild and severe cases, respectively (Sup- mum free energy (MFE) secondary structure of the plementary Table S1). spike gene (positive sense for A > G and negative Most haplotypes per amplicon (1130/1536, 73.6%) sense for T > C) was generated for each mutant were present at ≤0.5% in the respective quasispecies using a loop-based energy model and a dynamic pro- (Supplementary Table S2). In general, each haplotype gramming algorithm [37], and the RNA folding tools had a single mutation relative to the master haplotype offered by the University of Vienna website resources (Supplementary Table S3). [38, 39]. Analysis of the quasispecies diversity indices in relation to COVID-19 severity showed that diversity (rare haplotype load at 1%, Shannon entropy, and Characterization of fixed mutations worldwide mutation frequency) tended to be higher in severe To study the type of mutations fixed in SARS-CoV-2 than in mild cases (Table 2). This trend could be isolates worldwide, we extracted the main SNVs per related to the effects of severe disease, to the higher position from the 353,341 events observed in 48,625 viral load, or to the lengthier disease duration in severe GISAID genomes by Mercatelli and Giorgi [40]. cases. However, some patients with mild disease Mutations supported by at least 10 genomes were con- showed high diversity, whereas some with severe dis- sidered the main SNVs. This computation eliminated ease showed low diversity (Supplementary Tables S1– prevalence effects due to bias in the geographic ori- S4) with no correlations with viral load, indicating that gins, resulting from the very different number of gen- higher diversity may be driven by the duration of the omes from different places. The aim of this infection. computation was not to identify the most prevalent SNVs, but instead, to identify any SNV observed in any genome worldwide and to obtain a count of Spike gene variants mutation types. Comparison of spike consensus sequences from all patients with sequences from the Wuhan reference genome (MN908947.3) found that all patients but Results one (P02) had the A1841G substitution (D614G), which classifies the virus as clade G, and two other Spike gene quasispecies diversity patients (P07 and P17) had the G3707T (C1236F) The SARS-CoV-2 spike gene was extracted, amplified, and G3231A (syn.) substitution, respectively. and sequenced in all patients. In total, 41,326,097 reads passed our quality filter, with a median coverage Minority variants of 145,391 reads per amplicon (interquartile range [IQR] 115,576–201,833). Coverage of amplicons In the total of patients, variant frequency was lower from patients with severe disease tended to be higher than 2% in most patients with mild disease except than that of patients with a mild condition (Table 2). for C0890T (S0297L) in P14 and C1467T (Syn) in In total, 1536 haplotypes were detected overall with P15, present at 2.4% and 2.35%, respectively (Sup- no large differences between the two clinical groups: a plementary Tables S4 and S5). Variants accounting median of 4 and 3.5 haplotypes per amplicon was for >5% of the viral quasispecies were seen in three Emerging Microbes & Infections 1781 Figure 1. Single nucleotide variants (SNVs) according to their frequency and position in the spike (S) region. Top panel shows SNVs with frequencies above 1%. Bottom panel shows SNVs below 1%. The spike regions are depicted in different background colours. CP, cytoplasmic domain; FP, fusion peptide; HR, heptad repeat; NTD, N-terminal domain; QS, quasispecies; RBD, receptor-binding domain; SP, signal peptide; TM, transmembrane. patients with severe disease (Figure 1 and Supplemen- bases in the spike gene), with 1090 SNVs. On determi- tary Table S4). In two of the three, variants were present nation of their abundance in the viral quasispecies, 6 at >10%: T2875C (Syn) at 21.85% in P08, and G3306T positions had SNVs at ≥5% (Figure 1 and Supplemen- (W1102C) at 11.6% and G1915T (G639C) at 10.8% in tary Table S4). Of the 1090 SNVs identified, 942 P11. The majority of highly frequent mutants (>2% of (86.4%) comprised haplotypes present in <0.5% of quasispecies) were found in severe cases (P08, P09, the quasispecies. Among these, 1052 (96.5%) were P10, P11, P12, P13, P16, P17) except for the two transitions, mainly A → G and T → C (908, 83.3%) mutations mentioned above in P14 and P15. (Figure 2, Table 3). The SARS-CoV-2 spike gene is Aggregation of all variants found in all patients by AT-rich (62.7%), containing 29.41% A and 33.25% frequency bin (Table 3) and by substitution type T. However, 83.3% of all substitutions involved A → yielded 642 polymorphic sites (16.8 8% of the 3822 G and T → C changes (ADAR-type) showing a clear 1782 J. Gregori et al. Table 3. Counts of observed SNVs by frequency bins expressed in % ranges, among the 18 patients. Frequency bins SNV (0,0.2] (0.2,0.5] (0.5,1] (1,2] (2,5] (5,10] (10,50] Sum % A->C 2 1 1 4 0.37% A->G 160 243 20 10 7 1 441 40.46% A->T 1 1 1 3 0.28% C->A C->G C->T 15 51 17 12 3 1 99 9.08% G->A 8 26 7 1 2 1 45 4.13% G->C 2 2 4 0.37% G->T 2 5 2 5 2 2 18 1.65% T->A 3 3 0.28% T->C 169 255 35 5 2 1 467 42.84% T->G 4 1 1 6 0.55% Sum 355 592 83 38 16 3 3 1090 % 32.57% 54.31% 7.61% 3.49% 1.47% 0.28% 0.28% The frequencies are as observed in the corresponding amplicon quasispecies. Transitions (in bold) accounted for 96.44% of observed SNVs. ADAR-type substitutions, A → G and T → C, accounted for 84.07% of the SNVs in this population. Most SNVs (86.92%) show abundances below 0.5% in the corre- sponding quasispecies. bias with respect to random substitutions. Non- values per patient were, respectively, 1124 A, 723 C, synonymous substitutions dominated in all abun- 704 G, 1271 T, 37 polymorphic sites (0.97%), 10 A dance bins except (5%, 10%]. The ratio of transitions → G (0.89% of A), and 12 T → C (0.94% of T). Overall to transversions and the ratio of ADAR-type to non- there were 1,090 SNVs, 441 (40.5%) A → G tran- ADAR-type substitutions varied considerably as the sitions, and 467 (40.5%) T → C transitions. abundance of SNVs increased, with values of 172.0 On analysis of the number of substitutions in relation and 14.0, respectively, in the lowest frequency bin to their sites in the alignment of all 18 patients (Sup- (0%, 0.2%], and 0.5 and 0.5, respectively, in the highest plementary Table S9), we found 441 SNVs with A → bin (10%, 50%] (Table 4). Only 6 SNVs were present at G in any of the 230 sites with this type of substitution, a frequency of >5% (Figure 1, panel A, Supplementary and 467 SNVs with T → C in any of the 254 sites with Table S4). The percentage of most SNVs was very low this type of substitution. Conversely, there were 45 (<1%) (Figure 1, panel B). The fact that the ADAR- SNVs with G → A in any of the 41 sites with this type type to non-ADAR-type ratio was significantly higher of substitution, and 99 SNVs with T → C in any of in the lowest bins suggests that ADAR1 was the cause the 88 sites with this type of substitution. All the trans- of those substitutions, with hypermutation effects, versions observed were unique. In fact, most SNVs whereas the higher bins would incorporate the involved a single haplotype (Supplementary Table S3). effects of haplotype fitness and selection (Table 4, Sup- plementary Figures S1 and S2). Analysis of the mutations in relation to their pos- Concurrent SNVs ition in the codon showed that the A → G and T → Analysis of the quasispecies to detect SNVs common C transitions were randomly distributed and not to more than one patient found that some variants restricted to specific positions (Figure 2). A significant were present in several patients, including one overrepresentation of (A → G)/A relative to (G → observed in 11. The mutations common to various A)/G, and of (T → C)/T relative to (C → T)/C, both patients showed no predilection for any specific region overall and in each codon position, was confirmed by of the spike. The A2194G mutation was detected in 11/ Fisher tests (Supplementary Table S6). The homogen- 18 patients at a minimal abundance of 0.18% and a eity of (A → G)/A substitutions with respect to (T → maximum of 0.26% (Table 5). All concurrent C)/T could not be rejected by Fisher test, which is con- mutations were present at <1.0% (median [IQR] = sistent with the hypothesis that As are equally-likely 0.21% [0.19; 0.245]). Given the low abundance of edited on both strands (See Supplementary Table S7). these concurrent variants along with the very high Among the transversions, there was also a significant abundance of the dominant haplotypes (>90%), we overrepresentation of (G → T)/G relative to (T → G)/T. can assume that the former were generated indepen- The analysis was completed by computing the fol- dently in each patient during viral replication lowing values for each patient: the number of bases (Table 5). Moreover, all master sequences had abun- in the respective consensus sequence, the number of dances well above 50%, with most of them above polymorphic sites, the number of poly sites with A 90%; hence, amplicon master sequences were equival- → G transitions, the number of poly sites with T → ent to the amplicon quasispecies consensus sequences. C transitions, and the corresponding fractions with Because of the high abundance of dominant haplo- respect to the total of sites, the total of A sites, and types, we suggest that each patient was infected by a the total of T sites (Supplementary Table S8). Median single haplotype (bottlenecking effect), and that all Emerging Microbes & Infections 1783 Figure 2. Changes observed between pairs of bases by codon position. In abscissas, the WT base, with variants in the bars. In ordinates, the ratio (#X → Y substitutions) / (#X sites). The pattern of substitutions is homogeneous in all three codon sites, despite the A → G and T → C biases observed. The incidence of transversions was very low, with G → T and T → A dominating. the variants observed had been generated indepen- substitutions (A1988G and A3785G) took place in dently within each patient. the helix, one in the hairpin (A2194G), and one in a The level of concurrency may be better understood multibranched loop (A2324G) (Figure 4, panel A-D). if we consider that overall there were 441 A → G sub- In contrast, T → C changes occurred mainly in stitutions in different haplotypes from the 18 patients, internal hairpins, with the exception of T1389C, 436 when SNV duplicates in the same patient were which was found in the helix (Figure 4, panels E-H). removed, corresponding to 230 unique sites, with a None of the substitutions significantly changed the ratio of 436/230 = 1.9. Furthermore, there were 467 RNA secondary structure except for T1683C, which T → C substitutions in different haplotypes, 460 changed the internal loop into a predicted bulge (Sup- when duplicates in the same patient were removed, plementary Figure S3). corresponding to 254 unique sites, with a ratio of Inspection of the dinucleotides immediately before 460/254 = 1.8. and after the A → G substitution sites showed a promi- nent 18.8% of AG and 19.5% of TG at 5’, and 25.7% of CA at 3’ compared to the randomly expected frequency Substitution context in the spike gene of 6.25% (1/16). The same analysis on the T → C substi- secondary structure tution sites also showed prominent dinucleotides: 19.8% TG, 14.1% CC, and 12.1% GC at 5’, and 18.9% CA, In all 18 samples, A → G substitutions were studied to 17.4% TT, and 15.7% CT at 3’ (Supplementary Figure determine whether they were affected by neighbouring S4). This bias may be partially explained by the rep- nucleotides; no specific pattern was detected at either etition of some substitutions among the patients studied. the 5’ base (16.7% A, 6.0% T, 49.3% G, 28.0% C) or the 3’ base (6.4% A, 7.3% T, 36.0% G, 50.2% C). However, analysis of sites with recurrent mutations found in 5 or Fixed mutations worldwide more patients showed a more specific pattern at the 5’ base (5.6% C, 94.4% G) and 3’ base (69.4% C, 30.6% G) SNVs observed in 48,625 GISAID genomes worldwide (Figure 3). by Mercatelli and Giorgi (2020) [40] were identified. Analysis of the substitution context of the RNA sec- Among the SNV types, we found that the fixed ondary structure only in SNVs common to at least 8 mutations were mainly C → T and G → T rather than patients showed that the A → G transition in two ADAR1-associated changes (Supplementary Table S10). Table 4. Ratio of transitions to transversions, ADAR-like substitutions to other substitution types, and synonymous to non- synonymous substitutions as observed by frequency bins. Frequency bins (0,0.2] (0.2,0.5] (0.5,1] (1,2] (2,5] (5,10] (10,50] Tr/Tv 172 33.2 26.0 2.70 7.00 Inf 0.50 ADAR/other 14 5.54 2.12 0.61 1.29 0.50 0.50 Syn/NonSyn 0.70 0.52 0.65 0.68 0.33 2.00 0.50 These results are obtained from Table 3. 1784 J. Gregori et al. Table 5. Concurrent SNVs. Pos WT Var N Min Frq (%) Max Frq (%) 2194 A G 11 0.18 0.26 1389 T C 10 0.18 0.23 2324 A G 10 0.18 0.40 1683 T C 9 0.18 0.23 1988 A G 9 0.18 0.29 512 T C 8 0.18 0.27 2092 T C 8 0.20 0.28 3785 A G 8 0.18 0.43 1224 A G 7 0.19 0.26 1233 T C 7 0.19 0.24 1616 T C 7 0.18 0.32 2443 A G 7 0.18 0.54 2942 T C 7 0.22 0.55 467 A G 6 0.18 0.34 885 T C 6 0.18 0.27 1147 T C 6 0.20 0.26 2479 A G 6 0.20 0.26 2902 T C 6 0.19 0.24 2913 T C 6 0.18 0.24 2963 A G 6 0.18 0.29 3298 A G 6 0.19 0.24 3398 T C 6 0.18 0.34 3452 A G 6 0.18 0.28 3588 T C 6 0.18 0.31 3779 A G 6 0.19 0.30 Position, wild type base, variant base, number of patients with the substi- tution, and minimum and maximum frequency among all haplotypes showing the substitution. One SNV was seen in 11 patients, 2 were seen in 10 patients, 2 in 9 patients, 3 in 8 patients, 5 in 7 patients, 11 in 6 patients, as well as 13 in 5 patients and 25 in 4 patients that are not in the Table. Discussion Although SARS-CoV-2 RNA replication has a proof- reading control, deletions and mutations have been reported in the spike gene [10]. NGS enables the deep study of the viral quasispecies and identification of mutations that might be selected and affect viral replication or pathogenesis. In this study, we analyzed samples from a group of COVID-19 patients with mild or severe disease by NGS at the start of the pan- demic to investigate the complexity and variability of the viral quasispecies in the SARS-CoV-2 spike gene. The variant analysis showed that all overlapping amplicons in the 18 patients had a prominent, highly abundant haplotype as well as a few haplotypes at low abundance bearing a single variant in most cases. Despite the deep analysis performed, no haplo- type present in >0.8% of the quasispecies was found to be common to more than one patient. This finding supports the idea that only the most abundant haplo- types were transmitted person-to-person. Further- more, it would indicate that there is a low probability that variants would be generated along the transmission chain. Instead, variant dominance would mainly occur in patients with a lengthy period of viremia. The diversity analysis results suggest that SARS- CoV-2 quasispecies maturity and diversity are mainly driven by the duration of the infection, rather than by clinical severity. We found that SARS-CoV-2 is less Figure 3. Location of concurrent substitutions occurring in diverse than other RNA viruses, such as HCV, HIV, more than 8 patients. A–D shows substitutions A to G and E–H shows substitutions T to C in the predicted complete and HBV, and that it makes use of mechanisms of spike RNA gene secondary structure of the positive strand in the wild-type (WT) genome. Emerging Microbes & Infections 1785 Figure 4. Sequence context of SARS-CoV-2 RNA edited sites. Top left includes all sites with A → G transitions. Top right includes sites with A → G concurrent transitions occurring in 4 or more patients. Bottom left includes all sites with T → C transitions. Bot- tom right includes sites with T → C concurrent transitions occurring in 4 or more patients. variability other than its own RNA polymerase. The individuals. In this situation, immune escape or resist- large coronavirus genome requires a balance between ant mutants may arise and be selected. This may have stability and diversity. The proofreading exoribonu- been the case of the British (VOC 202012/01), Brazi- clease 3’-5’ activity by nsp14 might have an impact lian (P.1 from B.1.1.28), Japanese (B.1.1.248), and on the high fidelity of replication [1, 41, 42]. However, South African (501Y.V2) variants, all sharing the viruses can exploit other sources of variation such as N501Y mutation in the receptor-binding domain of recombination [4] and acceptance of deletions [4– the spike protein. Transmission of minority variants 10] and insertions. In addition, they can use cellular could, however, also occur in patients with lengthy mechanisms of innate immunity, such as ADAR1 edit- infection and a high viral load, where there would be ing activity [43]. Previous studies investigating the dis- time enough for the virus to replicate and generate a tribution of SARS-CoV-2 mutations [40, 44] are based highly fit variant, or for a variant to occur by chance on consensus sequences downloaded from the after a bottlenecking event [48]. However, in the pre- GISAID repository database [45, 46]; hence, viral gen- sent study viral load did not correlate with COVID-19 ome subpopulations within the quasispecies of indi- disease severity. Severity has been linked with the vidual patients cannot be analyzed. Here we found patients’ demographic data, such as age and sex, and that 86.9% of the variants observed using very deep particularly with pre-existing comorbidities such as sequencing (>70,000 reads) were present at <0.5%; hypertension, diabetes, cardiovascular disease, chronic hence, they can be considered rare variants. These lung disease, and cancer, which have been associated rare variants appeared at the onset of nucleotide edi- with greater severity and a higher fatality rate. Fur- tion, prior to any fitness or fixation effects. Rare viral thermore, COVID-19 contributes to cardiovascular haplotypes may be more representative of the under- complications, including acute myocardial injury as lying mechanism causing the substitutions, whereas a result of acute coronary syndrome, myocarditis, mutants observed at higher rates additionally rep- stress cardiomyopathy, arrhythmias, cardiogenic resent the selected haplotypes with the highest shock, and cardiac arrest [49, 50]. fitness [47]. Accordingly, rare variants will mainly be Most variants found (96.4%) were transitions, transient because of their low abundance, low fitness mainly A → G and T → C substitutions (88.1%), con- value, or limited time to replicate and reach a signifi- sistent with the action of ADAR1 [23]. ADARs edit cant load, unless a selective pressure (e.g. treatment specific sites or hyper-mutate full viral genes or gen- with convalescent sera or antiviral monotherapy) omes [51], which leads to mutations that can be occurs in patients with lengthy persistence of SARS- selected depending on the dynamic interaction CoV-2 infection, such as immunocompromised between the virus and host. Several viruses have 1786 J. Gregori et al. used this genetic mechanism to generate variability only a single substitution (Supplementary Table S4). that supports increases in infectivity and evolutionary Nonetheless, we should take into account that each potential [16–18, 43, 51, 52]. Various SNVs observed amplicon was roughly a 1/13th of the spike gene, in this study were common to more than one patient, which in turn is the most variable region, but only although at a low rate (<1%), suggesting that there are about 10% of the full genome. Hence, if ADAR intro- preferred nucleotide positions or “hot spots” within duced only 2 substitutions on edited spike genes, it the genome where host enzymes such as ADAR1 can would imply around 20 mutations per genome, most act. The probability that a particular point mutation of them non-synonymous (Table 4), which would would be found in multiple patients (hypothesis of suffice to impair the edited genomes. Regardless of independence) was almost negligible. Of note, all sub- the mechanism involved, the data suggest that stitutions present in >3 patients were either A → G or ADAR has an antiviral action on the spike gene, T → C (at nt31 and 32, respectively), suggesting that although to a limited extent, as most of the replicating ADAR1 had a role in generating them. It is important virions remained unaffected. to highlight that the number of true substitutions ADAR1 acts on dsRNA, and editing by this enzyme should be even higher than the values reported in could occur in either the positive or the negative this study because of the very stringent filters used strand. ADAR1 editing of the positive strand during to avoid sequencing artifacts. viral replication, in which dsRNA is formed, should To determine whether ADAR1 editing is affected lead to an A → I transition, whereas ADAR1 editing by the 5’ or 3’ neighbours of the substituted nucleotide of the negative strand would result in a T → C tran- [36], all A → G and T → C substitutions were ana- sition, once the negative strand has served as the tem- lyzed, but no pattern emerged indicating a preferred plate for a new positive-strand [55]. Because of its neighbour on either the 5’ or 3’ end base. However, considerable length, SARS-CoV-2 RNA has a complex when only concurrent substitutions common to four secondary structure, with many double strands form- or more patients were analyzed, the pattern was sur- ing helices, hairpins, bulges, and loops [37, 38, 56]. prisingly more specific on both end bases. The 5’ Study of the spike secondary RNA structure in substi- neighbour preference of mutated A → G was G fol- tutions common to more than 8 patients showed that lowed by C and A, whereas the 3’ neighbour prefer- most mutations occurred in hairpins or loops. This ence of mutated T → C was C followed by T and G. suggests that the editing preferably occurred during Neither the high A → G or T → C bias nor the replication. These results support the concept of mul- neighbour preference led to the fixed mutations in the tiple ADAR1 editions in a single genome, in either the consensus genomes of the GISAID [45, 46]. Most positive or the negative strand, during RNA replica- recurrent substitutions occurring in four or more of tion. Hypermutated genomes could carry deleterious our patients (51 of 63) could not be found in the mutations compromising fitness, and in a scenario GISAID genomes analyzed. Nine were found in a single of low selective pressure, these mutants would have GISAID genome, two were found in two GISAID gen- a very low probability of being selected and trans- omes, and only one in three GISAID genomes. In our mitted. However, ADAR1 editing is a source of vari- view, the most likely explanation for this would be ation providing a reservoir of genomes that have a related to the multiple non-specific editions (hypermu- chance to be selected in response to an environmental tations) produced by ADAR. A small fraction of repli- challenge (e.g. immune response, vaccine, antiviral cating genomes underwent ADAR edition and these treatment), especially in long-term infected patients. sequences were edited several times, reducing the like- This study has some limitations, in particular, the lihood of replication fitness. That is why most ADAR- small number of patients included, resulting from type substitutions were detected in genomes present at the need for deep sequencing of each sample, which low rates (genomes with poor fitness) and why very few generates an enormous amount of data. In addition, ADAR1-type substitutions were fixed and detected in tissue samples from the patients studied were not consensus genomes (Supplementary Table S10). Thus, available, and it was impossible to measure intracellu- a limited ratio of available ADAR molecules to replicat- lar ADAR-1 activity. However, the findings open the ing genomes (stoichiometric reason), and a longer time door to investigation of whether new full variants of of engagement of each ADAR molecule due to the large concern may arise because of the effects of ADAR1 genome size, would limit the extension of ADAR1 on the quasispecies. hypermutation activity to a larger number of genomes, In conclusion, SARS-CoV-2 quasispecies variability as has been observed in smaller viruses [53, 54]. Our was higher in patients with long-lasting infection, results suggest that most fixed mutations would come regardless of their severity status. Deep sequencing from replication errors due to the polymerase, despite showed that most SARS-CoV-2 mutations present in the proofreading mechanism. the quasispecies were of the type produced by the The ADAR hypermutation hypothesis could be action of ADAR1 editing enzymes, rather than substi- questioned by the fact that most haplotypes carried tutions inserted by the RNA-dependent RNA Emerging Microbes & Infections 1787 polymerase (nsp12). Furthermore, ADAR1 editing Josep Francesc Abril http://orcid.org/0000-0001-7793- preferably occurred during replication (dsRNA with 589X a+ and −strand), as the effects were similar on both Maria Gema Codina http://orcid.org/0000-0002-5313- 5086 strands. Nucleotide editing by ADAR1 took place on Ariadna Rando http://orcid.org/0000-0003-4555-7286 a small fraction of replicating genomes, but often Juliana Esperalba http://orcid.org/0000-0003-1326-1341 resulted in multiple editions in the same genome, Elena Sulleiro http://orcid.org/0000-0002-9783-6060 which can compromise genome fitness and replica- Joan Joseph http://orcid.org/0000-0002-4106-0208 tion. Further study at the long term showed very few Narcís Saubí http://orcid.org/0000-0001-5493-1044 Sergi Colomer-Castell http://orcid.org/0000-0002-0307- ADAR-like mutations fixed in the GISAID consensus 7455 genomes, which is consistent with the idea that Mari Carmen Martin http://orcid.org/0000-0003-1327- ADAR1 edition acts in the manner of a moderate anti- 8813 viral by causing hypermutations that impair viral Juan Ignacio Esteban http://orcid.org/0000-0001-5085- fitness. Finally, based on our data, we suggest that 917X most, if not all, SARS-CoV-2 mutations that can be Tomas Pumarola http://orcid.org/0000-0002-5171-7461 Francisco Rodriguez-Frias http://orcid.org/0000-0001- transmitted and become fixed in the consensus 9058-4641 sequence are likely nsp12-induced and are produced Andrés Antón http://orcid.org/0000-0002-1476-0815 in patients with high viremia and lengthy infection. Josep Quer http://orcid.org/0000-0003-0014-084X Acknowledgments The authors thank Celine Cavallo for English language sup- References port. JG and MFC performed the ADAR1 deep data analy- sis, the graphic analysis and composed the figures and [1] Eckerle LD, Becker MM, Halpin RA, et al. Infidelity of tables; MP, CCampos, CG-C, CA, SC-C, MCM and CCas- SARS-CoV Nsp14-exonuclease mutant virus replica- tillo conducted the experiments, performed the RNA extrac- tion is revealed by complete genome sequencing. tion, PCR and sequencing and the RNA secondary analysis; PLoS Pathog. 2010;6:e1000896. MGC, AR, JE, ES collected patients samples, patients data, [2] Peck KM, Lauring AS. Complexities of viral mutation and performed patient follow up; JJ and NS helped with rates. J Virol. 2018;92; doi:10.1128/JVI.01031-17. the theoretical analysis of data and the RNA secondary [3] Rambaut A, Holmes EC, O’Toole Á, et al. A dynamic analysis; FR-F, TP and JIE revised the manuscript; AA nomenclature proposal for SARS-CoV-2 lineages to and JQ directed the work, got the funding, analyzed the assist genomic epidemiology. Nat Microbiol. 2020. data, and wrote the manuscript with the help of JG. doi:10.1038/s41564-020-0770-5 [4] Rowe CL, Fleming JO, Nathan MJ, et al. Generation of coronavirus spike deletion variants by high-frequency Disclosure statement recombination at regions of predicted RNA secondary structure. J Virol. 1997;71:6183–6190. No potential conflict of interest was reported by the author [5] Rowe CL, Baker SC, Nathan MJ, et al. Evolution of (s). mouse hepatitis virus: detection and characterization of spike deletion variants during persistent infection. J Virol. 1997;71:2959–2969. Funding [6] Muth D, Corman VM, Roth H, et al. Attenuation of replication by a 29 nucleotide deletion in SARS-coro- This study was supported by the Direcció General de navirus acquired during the early stages of human-to- Recerca i Innovació en Salut (DGRIS) of the Catalan Health human transmission. Sci Rep. 2018;8:15177. Ministry, Generalitat de Catalunya, through the Vall d’Heb- [7] Lau S-Y, Wang P, Mok BW-Y, et al. Attenuated SARS- ron Institut de Recerca (VHIR); the European Regional CoV-2 variants with deletions at the S1/S2 junction. Development Fund (ERDF) “A way to achieve Europe” by Emerg Microbes Infect. 2020;9:837–842. the Spanish Network for Research in Infectious Diseases [8] Wong YC, Lau SY, Wang To KK, et al. Natural trans- [grant number REIPI RD16/0016/0003]; Instituto de Salud mission of bat-like SARS-CoV-2 deltaPRRA variants Carlos III [grant number FIS PI19/00301]; and Centro in COVID-19 patients. Clin Infect Dis. 2020. doi:10. para el Desarrollo Tecnologico Industrial (CDTI) from the 1093/cid/ciaa953 Ministry of Economic Affairs and Digital Transformation [9] Young BE, Fong S-W, Chan Y-H, et al. Effects of a major [grant number IDI-20200297]. deletion in the SARS-CoV-2 genome on the severity of infection and the inflammatory response: an observa- ORCID tional cohort study. Lancet. 2020;396:603–611. [10] Andrés C, Garcia-Cehic D, Gregori J, et al. Naturally Josep Gregori http://orcid.org/0000-0002-4253-8015 occurring SARS-CoV-2 gene deletions close to the Maria Francesca Cortese http://orcid.org/0000-0002- spike S1/S2 cleavage site in the viral quasispecies of 4318-532X COVID19 patients. Emerg Microbes Infect. Maria Piñana http://orcid.org/0000-0002-4766-2613 2020;9:1900–1911. Carolina Campos http://orcid.org/0000-0002-0132-7027 [11] Majzoub K, Wrensch F, Baumert TF. The innate anti- Damir Garcia-Cehic http://orcid.org/0000-0002-0009- viral response in animals: an evolutionary perspective 038X from flagellates to humans. Viruses. 2019;11; doi:10. Cristina Andrés http://orcid.org/0000-0002-3200-0895 3390/v11080758. 1788 J. Gregori et al. [12] Cattaneo R, Schmid A, Eschle D, et al. Biased hyper- [31] Morgan M, Anders S, Lawrence M, et al. Shortread: a mutation and other genetic changes in defective bioconductor package for input, quality assessment measles viruses in human brain infections. Cell. and exploration of high-throughput sequence data. 1988;55:255–265. Bioinformatics. 2009;25:2607–2608. [13] Lehmann KA, Bass BL. The importance of internal [32] Wagih O. Ggseqlogo: a versatile R package for drawing loops within RNA substrates of ADAR1. J Mol Biol. sequence logos. Bioinformatics. 2017;33:3645–3647. 1999;291:1–13. [33] Gallego I, Gregori J, Soria ME, et al. Resistance of high [14] Samuel CE. ADARs: viruses and innate immunity. fitness hepatitis C virus to lethal mutagenesis. Curr Top Microbiol Immunol. 2012;353:163–195. Virology. 2018;523:100–109. [15] Piontkivska H, Matos LF, Paul S, et al. Role of host-dri- [34] Gregori J, Perales C, Rodriguez-Frias F, et al. Viral ven mutagenesis in determining genome evolution of quasispecies complexity measures. Virology. sigma virus (DMelSV; Rhabdoviridae) in Drosophila 2016;493. doi:10.1016/j.virol.2016.03.017. melanogaster. Genome Biol Evol. 2016;8:2952–2963. [35] Gregori J, Soria ME, Gallego I, et al. Rare haplotype [16] Phuphuakrat A, Kraiwong R, Boonarkart C, et al. load as marker for lethal mutagenesis. PLoS One. Double-stranded RNA adenosine deaminases enhance 2018;13:e0204877. expression of human immunodeficiency virus type 1 [36] Lehmann KA, Bass BL. Double-stranded RNA adeno- proteins. J Virol. 2008;82:10864–10872. sine deaminases ADAR1 and ADAR2 have overlap- [17] Tang Y-D, Na L, Fu L-H, et al. Double-stranded RNA- ping specificities. Biochemistry. 2000;39:12875–12884. specific adenosine deaminase 1 (ADAR1) promotes [37] Zuker M, Stiegler P. Optimal computer folding of large EIAV replication and infectivity. Virology. RNA sequences using thermodynamics and auxiliary 2015;476:364–371. information. Nucleic Acids Res. 1981;9:133–148. [18] Doria M, Neri F, Gallo A, et al. Editing of HIV-1 RNA [38] Gruber AR, Lorenz R, Bernhart SH, et al. The Vienna by the double-stranded RNA deaminase ADAR1 RNA websuite. Nucleic Acids Res. 2008;36:W70–W74. stimulates viral infection. Nucleic Acids Res. [39] Lorenz R, Bernhart SH, Höner Zu Siederdissen C, et al. 2009;37:5848–5858. Vienna RNA package 2.0. Algorith Mol Biol. 2011;6:26. [19] Casey JL. Control of ADAR1 editing of hepatitis delta [40] Mercatelli D, Giorgi FM. Geographic and genomic virus RNAs. Curr Top Microbiol Immunol. distribution of SARS-CoV-2 mutations. Front 2012;353:123–143. Microbiol. 2020;11:1800. [20] Clerzius G, Gélinas J-F, Gatignol A. Multiple levels of [41] Smith EC, Denison MR. Coronaviruses as DNA PKR inhibition during HIV-1 replication. Rev Med wannabes: a new model for the regulation of RNA Virol. 2011;21:42–53. virus replication fidelity. PLoS Pathog. 2013;9:e1003760. [21] Wong SK, Lazinski DW. Replicating hepatitis delta [42] Gribble J, Pruijssers AJ, Agostini ML, et al. The coro- virus RNA is edited in the nucleus by the small form navirus proofreading exoribonuclease mediates exten- of ADAR1. Proc Natl Acad Sci U S A. sive viral recombination. bioRxiv. 2020. doi:10.1101/ 2002;99:15118–15123. 2020.04.23.057786 [22] Rangan R, Zheludev IN, Das R. RNA genome conser- [43] Tomaselli S, Galeano F, Locatelli F, et al. ADARs and the vation and secondary structure in SARS-CoV-2 and balance game between virus infection and innate immune SARS-related viruses. bioRxiv Prepr Serv Biol. 2020. cell response. Curr Issues Mol Biol. 2015;17:37–51. doi:10.1101/2020.03.27.012906 [44] van Dorp L, Acman M, Richard D, et al. Emergence of [23] Di Giorgio S, Martignano F, Torcia MG, et al. genomic diversity and recurrent mutations in SARS- Evidence for host-dependent RNA editing in the tran- CoV-2. Infect Genet Evol J Mol Epidemiol Evol scriptome of SARS-CoV-2. Sci Adv. 2020;6:eabb5813. Genet Infect Dis. 2020;83:104351. [24] Chen Q, Perales C, Soria ME, et al. Deep-sequencing [45] Shu Y, McCauley J. GISAID: global initiative on shar- reveals broad subtype-specific HCV resistance ing all influenza data – from vision to reality. Euro mutations associated with treatment failure. Surveill Bull Eur sur les Mal Transm Eur Commun Antiviral Res. 2020;174. doi:10.1016/j.antiviral.2019. Dis Bull. 2017;22. doi:10.2807/1560-7917.ES.2017.22. 104694. 13.30494. [25] Perales C, Chen Q, Soria ME, et al. Baseline hepatitis C [46] Elbe S, Buckland-Merrett G. Data, disease and diplo- virus resistance-associated substitutions present at fre- macy: GISAID’s innovative contribution to global quencies lower than 15% may be clinically significant. health. Glob Challeng. 2017;1:33–46. Infect Drug Resist. 2018;11:2207–2210. [47] García-Crespo C, Soria ME, Gallego I, et al. Dissimilar [26] Magoc T, Salzberg SL. FLASH: fast length adjustment conservation pattern in hepatitis C virus mutant spec- of short reads to improve genome assemblies. tra, consensus sequences, and data banks. J Clin Med. Bioinformatics. 2011;27:2957–2963. 2020;9. doi:10.3390/jcm9113450. [27] Gregori J, Salicru M, Domingo E, et al. Inference with [48] Novella IS, Quer J, Domingo E, et al. Exponential viral quasispecies diversity indices: clonal and NGS fitness gains of RNA virus populations are limited by approaches. Bioinformatics. 2014;30; doi:10.1093/ bottleneck effects. J Virol. 1999;73:1668–1671. bioinformatics/btt768. [49] Kang Y, Chen T, Mui D, et al. Cardiovascular manifes- [28] R Core Team. R: a language and environment for stat- tations and treatment considerations in COVID-19. istical computing. Vienna, Austria: R Foundation for Heart. 2020;106:1132–1141. Statistical Computing; 2020. [50] Mapping the human genetic architecture of COVID- [29] Pages H, Aboyoun P, Gentleman R, et al. Biostrings: 19. Nature. 2021. doi:10.1038/s41586-021-03767-x string objects representing biological sequences, and [51] Cattaneo R. Biased (A–>I) hypermutation of animal matching algorithms. R package 2.38.4. 2012. RNA virus genomes. Curr Opin Genet Dev. [30] Paradis E. Schliep K. ape 5.0: an environment for mod- 1994;4:895–900. ern phylogenetics and evolutionary analyses in R. [52] George CX, John L, Samuel CE. An RNA editor, ade- Bioinformatics. 2019;35:526–528. nosine deaminase acting on double-stranded RNA Emerging Microbes & Infections 1789 (ADAR1). J Interf Cytokine Res Off J Int Soc Interf hypermutagenesis reaction. Nucleic Acids Res. Cytokine Res. 2014;34:437–446. 1995;23:2573–2578. [53] Martinez MA, Vartanian JP, Wain-Hobson S. [55] Carpenter JA, Keegan LP, Wilfert L, et al. Evidence for Hypermutagenesis of RNA using human immunodefi- ADAR-induced hypermutation of the Drosophila ciency virus type 1 reverse transcriptase and biased sigma virus (Rhabdoviridae). BMC Genet. 2009;10:75. dNTP concentrations. Proc Natl Acad Sci USA. [56] Vandelli A, Monti M, Milanetti E, et al. Structural 1994;91:11787–11791. analysis of SARS-CoV-2 genome and predictions of [54] Martínez MA, Sala M, Vartanian JP, et al. Reverse the human interactome. Nucleic Acids Res. transcriptase and substrate dependence of the RNA 2020;48:11270–11283.