key: cord-0728434-ydguoz5m authors: Koçhan, Necla; Eskier, Doğa; Suner, Aslı; Karakülah, Gökhan; Oktay, Yavuz title: Different selection dynamics of S and RdRp between SARS-CoV-2 genomes with and without the dominant mutations date: 2021-03-03 journal: Infect Genet Evol DOI: 10.1016/j.meegid.2021.104796 sha: b8969e187c6e7a0d6c3299a6771252c44e8af182 doc_id: 728434 cord_uid: ydguoz5m SARS-CoV-2 is a betacoronavirus responsible for the COVID-19 pandemic that has affected millions of people worldwide. Pharmaceutical research against COVID-19 and the most frequently used tests for SARS-CoV-2 both depend on the genomic and peptide sequences of the virus for their robustness. Therefore, understanding the mutation rates and content of the virus is critical. Two key proteins for SARS-CoV-2 infection and replication are the S protein, responsible for viral entry into the cells, and RdRp, the RNA polymerase responsible for replicating the viral genome. Due to their roles in the viral cycle, these proteins are crucial for the fitness and infectiousness of the virus. Our previous findings had shown that the two most frequently observed mutations in the SARS-CoV-2 genome, 14408C > T in the RdRp coding region, and 23403A > G in the S gene, are correlated with higher mutation density over time. In this study, we further detail the selection dynamics and the mutation rates of SARS-CoV-2 genes, comparing them between isolates carrying both mutations, and isolates carrying neither. We find that the S gene and the RdRp coding region show the highest variance between the genotypes, and their selection dynamics contrast each other over time. The S gene displays higher tolerance for positive selection in mutant isolates early during the appearance of the double mutant genotype, and undergoes increasing negative selection over time, whereas the RdRp region in the mutant isolates shows strong negative selection throughout the pandemic. COVID-19 is a pandemic caused by the SARS-CoV-2 betacoronavirus, of zoonotic origin, with human-to-human transmission capacity, and since the first observed cases in the Wuhan province of China (Chan et al., 2020; Riou & Althaus, 2020) , it has infected 45,968,799, with over 1 million recorded deaths (as of 6 November 2020), with 69% of cases and 78% of deaths in Europe and the Americas. While its primary characteristics are respiratory system distress and long-onset fever, long term symptoms include neuroinvasion (Li, Bai & Hashikawa, 2020; , cardiovascular complications (Kochi et al., 2020; Zhu et al., 2020) , and gastrointestinal and liver damage (Lee, Huo & Huang, 2020; Xu et al., 2020) . Alongside these secondary symptoms, the disease is also capable of asymptomatic community-wide transmission due to its high transmissibility (Wong et al., 2020) , creating a unique global health risk of urgent interest. As a result, the high amount of frequently updated data on viral genomes on databases such as GISAID (Elbe & Buckland-Merrett, 2017) and NextStrain (Hadfield et al., 2018) provides researchers with invaluable resources to track the evolution of the virus. SARS-CoV-2 has a linear, single-stranded genome, and encodes its own genomic replication proteins, instead of relying on host proteins. The replication proteins are cleaved from the Orf1a 2020). Owing to their role in maintaining replication fidelity and directly affecting the mutationselection equilibrium of RNA viruses, these proteins are key targets of study in understanding the mutation accumulation and adaptive evolution of the virus (Eckerle et al., 2010; Peng et al., 2020) . One of the key mutations in the SARS-CoV-2 RdRp coding region is the 14408C>T mutation, responsible for the P323L aminoacid substitution, frequently co-occurring with the 23403A>G mutation of the S gene, responsible for the D614G substitution. In our previous study, we examined the top 10 most frequent mutations in the SARS-CoV-2 RdRp, and identified that four of them are associated with an increase in mutation density in two genes, the membrane glycoprotein (M) and the envelope glycoprotein (E) which are under relatively low selective pressure, and mutations in these genes can be used to infer reduced replication fidelity (Eskier et al., 2020a) . The mutation that most strongly affects the mutation density was the 14408C>T mutation, predicating it as a genotype of interest. In our follow-up study, we also identified the trend of a positive correlation between the presence of the 14408C>T mutation and increased mutation density in regions under low selective pressure persisted beyond our initial study (Eskier et al., 2020c) . 2020 July 2019), with the parameters "-i FASTA -o PHYLIP input.file > output.file" (Hubisz, Pollard & Siepel, 2011) . Following all the filters and steps described above, a total of 34,393 genomes were used in the study (18,125 for UK, 16,268 for the US). We further divided the remaining isolate sequences into mutant (MT) and wildtype (WT) genotypes, with isolates carrying the 14408C>T and 23403A>G mutations classified as MT, while isolates carrying neither mutation classified as WT. Modeling the daily mutation rate J o u r n a l P r e -p r o o f In this study, we considered the conservative constant mutation rate model (Vega et al., 2004) , given as follows: where, the number of mutations (d) occurred in an isolate from its ancestor is calculated as the multiplication of the daily mutation rate (µ) and the time difference (t) between the isolate and the NC_045512.2 reference genome. Based on the reference genome, we calculated the number of mutations in each isolate and the time difference between the isolate and the NC_045512.2 reference genome. We then estimated the daily mutation rates using least square fitting and measured the goodness of fit using the adjusted R-square statistic. The rate of accumulation of synonymous mutations per synonymous site (dS), the rate of accumulation of non-synonymous mutations per non-synonymous site (dN) and their ratio (dN/dS), which is also known as w or / were calculated using PAML package (version 4.9j, available at http://abacus.gene.ucl.ac.uk/software/paml.html) (Yang, 2007) , with the parameters "runmode=-2, Codonfreq=2" in codeml.ctl (pairwise sequence comparison). In order to better understand the impact of the mutant genotype (hereafter referred to as MT) with 14408C>T and 23403A>G mutations on SARS-CoV-2 genome evolution and mutation rates, we used standard mathematical evolutionary biology models (Nei, Kumar & Nei, 2000; Vega et al., 2004) . Our results show that evolutionary forces affect mutant genomes differently compared to genomes without these two mutations. Specifically, the S gene and the RdRp coding region seem to be the most differentially affected, albeit in opposite directions. In this study, we focused on the mutational dynamics of SARS-CoV-2 in UK and US populations, two regions with the highest number of isolate sequences submitted to GISAID. In particular, as our previous findings pointed to differences in mutation profiles between the two groups, we compared isolates that carried the 14408C>T and the 23403A>G mutations to isolates that did not. Isolates that carried both mutations were categorized as "mutant" (MT) while the isolates that carried neither mutation were categorized as "wild-type" (WT). Isolates Journal Pre-proof that carried one mutation but not the other were removed from the analysis. We therefore divided the isolates into the following categories: UK-WT (3,412 genomes), UK-MT (14,713 genomes), US-WT (2,839 genomes) and US-MT (13,428 genomes). In order to identify how the mutation rates of the SARS-CoV-2 genes vary, we calculated the synonymous and non-synonymous mutation rates (dS and dN, respectively) for each gene and each country (UK and US) separately ( Figure 1 ). Our results show that the non-synonymous mutation rates are lower than the synonymous mutation rates for all genes except E. This indicates that most genes have been more strongly affected by negative selection than positive selection. To assess the deviations from neutral selection, we calculated the dN/dS ratio per day for each gene ( Figure 2 ). The dN/dS ratios for E, M, Nsp14, ORF6, ORF7a in MT and WT isolates do not display high disparity in either country. The two genomic regions with the highest disparity were the RdRp coding region and the S genes, two regions carrying the MT-genotype defining mutations. Moreover, while the dN/dS ratio decreased for the S gene in MT isolates in both countries, it increased for the RdRp gene in the same isolates, whereas the dN/dS ratio for the ORF7b and ORF8 genes show a significant increase in US-MT isolates compared to US-WT isolates, while a similar increase is not observed in the UK isolates. We also performed the Wilcoxon rank-sum test to investigate the difference between the genotypes for individual genes in each country (Table 1) . The results show a significant difference in dN, dS and the dN/dS ratio between WT and MT groups for both countries in the RdRp coding region and the S gene (p<0.001). These results also corroborate previous findings on S and RdRp, and underline the importance of these two proteins for SARS-CoV-2 evolution and fitness. Journal Pre-proof Our previous findings had shown a shift in terms of mutation density increase between two time periods (days 60-100 and days 101-140) in the UK and the US, as day 100 (2 April To validate the results obtained from the US isolates, we repeated our analyses with the UK isolate sequences. As the UK only had sufficient isolate sequences available for days 80-99 and 100-143, we focused on these two time periods. We calculated daily mutation rates, dN and dS rates, as well as the dN/dS ratios for each gene and genotype. Our results show a significant difference in dN and dS values and the dN/dS ratios of RdRp and S genes between WT and MT isolates in different time intervals (p<0.001) (Figure 4 ). For both genotypes, the dN and dS rates decreased significantly in the plateau stage for the RdRp region and the S gene. In addition, the dN/dS ratio for RdRp increased in days 100-143 in both WT and MT genotypes, whereas the dN/dS ratio of the S gene increased in MT isolates but decreased in WT isolates in days 100-143. In addition to the gene-level mutation rate analysis, we next sought to estimate the synonymous and non-synonymous mutation rates per site per day, and the daily dN/dS ratio for each gene and country, using the conservative model previously described by Vega (Berrio, Gartner & Wray, 2020) . In both the UK and the US, the dN/dS ratio is highly similar for the majority of genes across genotypes. MT isolates tend to have lower dN/dS ratios for most genes compared to WT isolates. It is noteworthy that, similar to the mutation analysis described above and in Figure 1 WT isolates compared to MTs, whereas the opposite is true for the S gene, with its dN/dS ratio being lower in WT isolates compared to MT isolates. It should be noted that Orf8 is the third most disparate gene between MT and WT isolates. The significance of this finding remains to be explored, as there is limited knowledge on the biological function of Orf8 and its possible relationship to the two mutations that define the MT genotype is not clear. However, its proposed role in immune evasion could possibly provide a plausible explanation. As the RdRp region and the S gene were the two genomic regions with the most different dN/dS ratios between MT and WT isolates, we modeled their mutation rates per site in the US for the four time periods described in Figure Spike mutation D614G has been under spotlight since May 2020 and now it is widely accepted as causing a more transmissible form of SARS-CoV-2. Less is known about its co-mutation in RdRp (P323L), however, it could be leading to altered mutational fidelity. We previously showed that the mutant viruses accumulate mutations differently compared to those without. In this study, we show that viruses with the dominant 23403A>G (S D614G) and 14408C>T (RdRp P323L) mutations also have lower dN/dS ratios compared to those without these two mutations, particularly at RdRp coding region and Orf8 gene. On the contrary, S gene, unlike the other genes, has higher dN/dS ratios in the mutant genomes. Moreover, temporal analyses point to opposite trends of selection for these two critical genes. While the S gene seems to be under stronger negative selection in WT genomes during the early stages of the pandemic, it is almost at equal levels between MT and WT genomes in the later stages. On the other hand, RdRp is under stronger overall negative selection in the MT genomes, particularly during the early stages. One could speculate that RdRp with the P323L (14408C>T) mutation in the MT genomes has been relatively fit, and did not undergo significant changes in J o u r n a l P r e -p r o o f its dN/dS ratio over time. On the other hand, as the pandemic progressed, S with the dominant D614G (23403A>G) mutation has become increasingly less tolerant to non-synonymous mutations. As there are no other frequently observed mutations in the S gene, this dramatic change can be explained by the presence of more limiting epidemiological conditions for the virus in the later stages that led to stronger negative selection over time. It should be noted that these two genes harbor two of the most frequently observed mutations in the SARS-CoV-2 genome, and it is likely that non-synonymous mutations in these two genes/gene regions have relatively greater impact on viral fitness. RdRp is one of the SARS-CoV-2 proteins responsible for replicating the viral genome, and it is expected that mutations in its sequence and structure might lead to a change in overall mutation rates, as we have shown in previous studies (Eskier et al., 2020b,c) . It is less clear how the D614G mutation in the S gene could affect the mutation rate; however, increased transmission across the population might lead to more permissive circumstances for novel mutations to arise. As there are very few sequenced isolates carrying the D614G mutation but not the P323L mutation, or vice versa, there is not yet sufficient information to assess the individual impact of the mutations on genome replication fidelity or viral fitness. As novel variants of SARS-CoV-2 that are more transmissible and resilient to extant vaccines and treatments arise, it is essential to discover and understand the mechanisms underlying the ongoing evolution of the virus. The double mutant D614G and P323L genotype significantly affects the evolution process, and mechanistic study of these mutations may lead to a clearer image of how and where future mutations might arise. Furthermore, our methodology and results highlight the importance of timely and widespread sequencing of viral isolates in the population, J o u r n a l P r e -p r o o f Journal Pre-proof as well as in-depth analysis of these sequences, to better predict the course of the pandemic. One question that remains to be answered is whether the RdRp 14408C>T mutation has any effect on differential mutation rates observed between the S gene and the RdRp coding region of the two genotypes. Finally, the underlying reason for the observed association between the dN/dS ratio of the Orf8 gene and the MT genotype remains to be investigated. Overall, our results point to S and RdRp as critical players of SARS-CoV-2 evolution. A better understanding of the molecular basis of the evolutionary differences between these two genes/gene regions and the rest of the genome, as well as the effects of S D614G and RdRp P323L mutations, could lead to more effective strategies in combating COVID-19. Declarations of interest: none J o u r n a l P r e -p r o o f Positive selection within the genomes of SARS-CoV-2 and other Coronaviruses independent of impact on protein function A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster Infidelity of SARS-CoV Nsp14-exonuclease mutant virus replication is revealed by complete genome sequencing Data, disease and diplomacy: GISAID's innovative contribution to global health RdRp mutations are associated with SARS-CoV-2 genome evolution Mutation density changes in SARS-CoV-2 are related to the pandemic stage but to a lesser extent in the dominant strain with mutations in spike and RdRp Mutations of SARS-CoV-2 nsp14 exhibit strong association with increased genome-wide mutation load Nextstrain: real-time tracking of pathogen evolution PHAST and RPHAST: phylogenetic analysis with space/time models MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform Structure of the SARS-CoV nsp12 polymerase bound to nsp7 and nsp8 co-factors Cardiac and arrhythmic complications in patients with COVID-19 Gastrointestinal and liver manifestations in patients with COVID-19 The neuroinvasive potential of SARS-CoV2 may play a role in the respiratory failure of COVID-19 patients Structural basis and functional analysis of the SARS coronavirus nsp14-nsp10 complex. the National Academy of Sciences of the Molecular Evolution and Phylogenetics The Curious Case of the Nidovirus Exoribonuclease: Its Role in RNA Synthesis and Replication Fidelity Structural and Biochemical Characterization of the nsp12-nsp7-nsp8 Core Polymerase Complex from SARS-CoV-2 Pattern of early human-to-human transmission of Wuhan A Structural View of SARS-CoV-2 RNA Replication Machinery: RNA Synthesis, Proofreading and Final Capping Mutational dynamics of the SARS coronavirus in cell culture and human populations isolated in 2003 Asymptomatic transmission of SARS-CoV-2 and implications for mass gatherings. Influenza and Other Respiratory Viruses Nervous system involvement after infection with COVID-19 and other coronaviruses Liver injury during highly pathogenic human coronavirus infections PAML 4: phylogenetic analysis by maximum likelihood Cardiovascular Complications in Patients with COVID-19: Consequences of Viral Toxicities and Host Immune Response Writing -Original Draft Writing -Original Draft; Gökhan Karakülah: Conceptualization, Methodology, Writing -Original Draft, Writing -Review & Editing, Supervision, Resources; Yavuz Oktay: Conceptualization, Methodology, Writing -Original Draft The authors would like to thank Mr. Alirıza Arıbaş from Izmir Biomedicine and Genome Center for his technical assistance. The authors also would like to extend their thanks to Izmir Biomedicine and Genome Center (IBG) COVID19 platform IBG-COVID19 for their support in implementing the study. The following grant information was disclosed by the authors:The Scientific and Technological Research Council of Turkey (TUBITAK-STAR) and Turkish Academy of Sciences Young Investigator Program (TÜBA-GEBİP).