key: cord-1007139-u14c6huo authors: Mourier, T.; Shuaib, M.; Hala, S.; Mfarrej, S.; Alofi, F.; Naeem, R.; Alsomali, A.; Jorgensen, D.; Subudhi, A. K.; Ben Rached, F.; Guan, Q.; Salunke, R.; Ooi, A.; Esau, L.; Douvropoulou, O.; Nugmanova, R.; Perumal, S.; Zhang, H.; Rajan, I.; Al-Omari, A.; Salih, S.; Shamsan, A.; Al Mutair, A.; Taha, J.; Alahmadi, A.; Khotani, N.; Alhamss, A.; Mahmoud, A.; Alquthami, K.; Dageeg, A.; Khogeer, A.; Hashem, A. M.; Moraga, P.; Volz, E.; Almontashiri, N.; Pain, A. title: Saudi Arabian SARS-CoV-2 genomes implicate a mutant Nucleocapsid protein in modulating host interactions and increased viral load in COVID-19 patients date: 2021-05-10 journal: nan DOI: 10.1101/2021.05.06.21256706 sha: 0a6f6e7d49eb681c3b44ae3f597e89b4c7f174c5 doc_id: 1007139 cord_uid: u14c6huo Monitoring SARS-CoV-2 spread and evolution through genome sequencing is essential in handling the COVID-19 pandemic. The availability of patient hospital records is crucial for linking the genomic sequence information to virus function during the course of infections. Here, we sequenced 892 SARS-CoV-2 genomes collected from patients in Saudi Arabia from March to August 2020. From the assembled sequences, we estimate the SARS-CoV-2 effective population size and infection rate and outline the epidemiological dynamics of import and transmission events during this period in Saudi Arabia. We show that two consecutive mutations (R203K/G204R) in the SARS-CoV-2 nucleocapsid (N) protein are associated with higher viral loads in COVID-19 patients. Our comparative biochemical analysis reveals that the mutant N protein displays enhanced viral RNA binding and differential interaction with key host proteins. We found hyper-phosphorylation of the adjacent serine site (S206) in the mutant N protein by mass-spectrometry analysis. Furthermore, analysis of the host cell transcriptome suggests that the mutant N protein results in dysregulated interferon response genes. We provide crucial information in linking the R203K/G204R mutations in the N protein as a major modulator of host-virus interactions and increased viral load and underline the potential of the nucleocapsid protein as a drug target during infection. the mutant N protein results in dysregulated interferon response genes. We provide crucial information in linking the R203K/G204R mutations in the N protein as a major modulator of host-virus interactions and increased viral load and underline the potential of the nucleocapsid protein as a drug target during infection. The emergence of novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes the respiratory coronavirus infectious disease 2019 (COVID- 19) , resulted in a pandemic that has triggered an unparalleled public health emergency 1,2 . The global spread of SARS-CoV-2 depended fundamentally on human mobility patterns. This is highly pertinent to a country like the Kingdom of Saudi Arabia, which as of 22nd February 2021 had a total of 374,691 cases and 6,457 deaths 3 . The kingdom frequently experiences major population movements, particularly religious mass gatherings. For instance, during Umrah and Hajj roughly 9.5 million pilgrims visit two Islamic sites in Makkah and Madinah annually 4, 5 and the Ministry of Health takes extraordinary public health measures to keep the pilgrims safe and major outbreaks have been by and large avoided in recent years. Further, an estimated 5 million Shiite Saudi nationals travel to Iran for pilgrimage, which became an early source of COVID-19 infections in the region 5, 6 . This movement has been reflected in the early phase of COVID-19 transmission within Saudi, as the first case was officially reported in Qatif (Eastern Region) on March 2 nd , 2020 7 . Genomic epidemiology of emerging viruses has proven to be a useful tool for outbreak investigation and tracking the pathogen's progress 8, 9 . Currently, over half a million complete and high coverage genomes are accessible on GISAID 10, 11 , which aids immensely in tracking the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint viral sequences globally 12 . Novel SARS-CoV-2 variants are continuously arising and besides providing signals for epidemiological tracking, a subset of the resulting variants will have a functional impact on transmission and infection [13] [14] [15] . It is therefore critical to monitor the genetic viral diversity throughout the pandemic. In this study, we sequenced 892 SARS-CoV-2 genomes from nasopharyngeal swab samples of patients from the four main cities, Jeddah, Makkah, Madinah, and Riyadh, as well as a small number of patients from the Eastern region of Saudi Arabia (Figure 1 , Table S1 , Table S2 ). We analyzed the genomes to investigate the nucleotide changes and multiple mutation events that represent the first 6 months of the locally circulating pandemic lineages of the SARS-CoV-2 in Saudi Arabia and searched for potential association of polymorphic sites in the genome with available hospital records including severe disease and case fatality rates among the COVID-19 patients. We performed phylogenetic analysis to visualize the genetic diversity of SARS-CoV-2 and the nature of transmission lineages during March-August, 2020. We have presented a snapshot of the genomic variation landscapes of the SARS-CoV-2 lineages in our study population and linked specific set of mutation events in the N gene to viral loads in a diverse population of COVID-19 patients in Saudi Arabia ( Figure S1 ). Finally, we experimentally show the functional impact of these mutations in the N protein on the virus' interactions with the host. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint (top) and the estimate of effective population size [Ne] , the relative population size required to produce the diversity seen in the sample (bottom) . The red horizontal red line represents an R of 1, the level required to sustain epidemic growth. Grey confidence areas denote the 95% credible intervals. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint We sequenced and assembled SARS-CoV-2 genomes from 892 patient samples. This group includes 144 patients that were placed in quarantine and had either mild symptoms or were asymptomatic. The remaining patients were all hospitalized (Table S2 ). Data on comorbidities were available for 689 patients with diabetes (39%) and hypertension (35%) being the most abundant (Table S3 ). Patient outcome data was available for 850 samples, and 199 patients (23%) died during hospitalization (Table S2) . From the 892 assembled viral genomes collected over a period of 6 months, we found a total of 836 single-nucleotide polymorphisms (SNPs) compared to the Wuhan SARS-CoV-2 reference (GenBank accession: NC_045512) ( Figure S2 ). The observed numbers of SNPs relative to the Wuhan reference follow the numbers observed in global samples ( Figure S3 ). We further detected 41 indels of which 26 reside in coding regions (Table S4 ). Most indels were specific to a single sample, and no identical indel was found in more than four samples. Compared with global SNP data, seven SNPs were found in higher frequencies (absolute difference > 0.1) in samples from Saudi Arabia ( Figure S2 ). These include the Spike protein D614G (A23403G) and three consecutive SNPs causing the R203K and G204R changes in the nucleoprotein (G28881A, G28882A, and G28883C). Together with all sequences from Saudi Arabia available on GISAID on December 31 st 2020, the assembled sequences were used to construct the effective population size and growth rate estimates of SARS-CoV2 over the course of the first wave of the epidemic. The skygrowth model 16 ( Figure 1D ) shows a downward trend in the effective reproduction number (R) over time with the timely introduction and maintenance of effective non-. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint pharmaceutical interventions by the Saudi Ministry of Health. Following the lifting of restrictions towards the end of June, the model estimates that R remained below or at 1 to the end of the period covered by the genetic data presented in this study. The effective population size (N e ) represents the relative diversity of the sequences collected in Saudi Arabia over the course of the outbreak ( Figure 1D ). The model predicts a peak in viral diversity at the beginning of June. This is ahead of the peak number of cases reported nationally and is likely influenced by the earlier peak in reported cases in the three cities, which contribute the most viral sequences to this analysis (Madinah, Makkah and Jeddah). A maximum-likelihood phylogenetic analysis revealed that samples from Saudi Arabia represent 5 major Nextstrain clades 12 , 19A-B and 20A-C ( Figure 2A ). This highlighted the clade 20A that all carried the Nucleocapsid (N) protein R203K/G204R mutations 17 with high incidences of ICU hospitalizations. These samples were predominantly coming from Jeddah. Through time-scaled phylogenies dates of importation events were then estimated for each clade. The majority of importations for all clades were inferred to have occurred early in the outbreak, primarily in March and early April ( Figure 2B ). Inferring importation events from a phylogenetic tree with estimated dating of nodes we see an early import from Asia followed by multiple imports from different continents ( Figure 2C ). A dated phylogeny of global samples showed that samples with the R203K/G204R SNPs are predominantly found in Nextstrain clades 20A, 20B and 20C, and do not form a monophyletic group ( Figure S4 ). Furthermore, a few samples are further found in the early appearing 19A and 19B clades. However, due to the limited number of mutations separating SARS-CoV-2 genomes . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint constructing a reliable and robust phylogeny is problematic 18 , and while different clades may be well supported, the exact relationship between clades is often less easily resolved. Although phylogenetic trees of SARS-CoV-2 genomes may appear to robustly reflect transmission events, collapsing branches with low support will typically result in extensive polytomies 19, 20 . Additionally, the placement of individual virus genomes may be hampered by systematic errors, homoplasies, potential recombination, or co-infection of multiple virus strains 18, 19, [21] [22] [23] . It is therefore not clear if the phylogenetic distribution of samples with R203K/G204R SNPs reflects multiple independent origins of the SNPs, although it is evident that the R203K/G204R SNPs appeared early in the pandemic spread ( Figure S4 ). Consistent with this, we find the earliest estimated importation events of R203K/G204R SNPs in late January 2020, most likely from Italy ( Figure S5 ). This thus suggest a slightly earlier importation date than the estimate of importation . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. South America . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 10, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint A genome-wide association study between SARS-CoV-2 SNPs and patient mortality identified the three consecutive SNPs (G28881A, G28882A, G28883C) underlying the R203K/G204R mutations ( Figure 3B , Figure S6 ). Of the 892 assembled genomes, 882 (98.9%) genomes either have the three reference alleles, GGG, or the three mutant alleles, AAC, at positions 28,881-28,883. This is similarly found in global samples deposited in GIASID in 2020, where 99.7% of samples with SNPs at positions 28,881-28,883 contain all three SNPs ( Figure S7 ). In our samples, no other SNPs co-occur with the R203K/G204R SNPs ( Figure S8 ). The frequency of the R203K/G204R SNPs is markedly higher in samples from Jeddah, where the observed frequency of 0.38 is more than 10-fold higher than the average of the other cities (Table S2) . Within-host polymorphism has been observed for the R203K/G204R SNPs either resulting from co-infection of multiple strains or cross-sample contamination 21 . Co-infection of SARS-CoV-2 is demonstrated through observations of recombination between genetically distinct lineages 26 . To rule out cross-sample contamination, we investigated the levels of within-host polymorphisms in a range of SNP positions and found this more consistent with cases of co-infection among patients rather than contamination issues (Supplementary Information, Figure S9 ). Using multivariable regression, we next evaluated the effect of the R203K/G204R SNPs on mortality, severity, and viral load in our COVID-19 patients samples for which limited amount of clinical meta-datasets were available. Disease severity was defined as deceased patients and patients admitted to ICU. For mortality and severity, we first fitted a linear model using R203K/G204R SNPs as a covariate. Then we fitted adjusted models by including gender, age, comorbidities, hospital and time. Additionally, the Spike protein D614G SNP that is associated . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint with higher viral load 13 was included. Age and time were included using smoothing splines to allow for potential non-linear relationships 27 . Using an unadjusted logistic regression, we observed a positive and statistically significant association between R203K/G204R SNPs and severity. Specifically, we found that the log-odds of severity increased by 1.16, 95% CI 0.70-1. 64 . In the adjusted model, the log-odds decreased to 0.66, 95% CI 0.01-1.32. That is, we found a borderline significant association between R203K/G204R SNPs and severity. The relationship between mortality and R203K/G204R SNPs was positive and statistically significant in the unadjusted model with log-odds equal to 1.39, 95% 0.99-1.79. We also found a positive and statistically significant relationship adjusting for D614G SNP, gender, age, comorbidities, and hospital (log-odds = 0.62, 95% 0.13-1.10). However, after adjusting for time as a variable, there was no longer any association between R203K/G204R SNPs and mortality (log-odds: 0.25, 95% CI -0.30-0.80). The models thus suggest a temporal component in our observations, and it is important to note that the recorded mortalities from Jeddah are concentrated on just a few dates ( Figure S10 ). Unfortunately, our data set does not allow us to assess if the observed mortality rates are the result of shifts in treatment regimes or admission procedures during the sampling window. We then tested if R203K/G204R SNPs were associated with higher viral copy numbers as indicated by the cycle threshold (Ct) values obtained through quantitative PCRs (see Methods). From the unadjusted regression we found a positive and statistically significant relationship between R203K/G204R SNPs and log 10 (viral copy number), with the mean of log 10 (viral copy number) values increasing by 1.03 units (95% CI 0.67-1.46). The significance was still observed in the adjusted model, although the relationship decreased to 0.57 units (95% CI 0.13-1.01) ( Figure 3C ). Similarly, the adjusted model showed a significant relationship between D614G . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint SNPs and log 10 (viral copy number), with the mean values increasing by 0.32 units (95% CI 0.13-1.01), consistent with earlier reports 13, 28 . The positive and statistically significant association of R203K/G204R SNPs with higher viral load in critical COVID-19 patients suggests its functional implications during viral infection. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint 29 The SARS-CoV-2 N protein binds the viral RNA genome and is central to viral replication 30 . Protein structure predictions have shown that the R203K/G204R mutations result in significant changes in protein structure 24 , theoretically destabilizing the N structure 31 , and potentially enhancing the protein's ability to bind RNA and alter its response to serine phosphorylation events 32 . The R203K/G204R mutations in the SARS-CoV-2 N protein are within the linkage region (LKR) containing the serine/arginine-rich motif (SR-rich motif) ( Figure 4A ), known to be involved in the oligomerization of N proteins 33, 34 . Protein cross-linking shows that N mutant protein (with the R203K/G204R mutations) has higher oligomerization potential compared to the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint control N protein (without the changed amino acids) at low protein concentration ( Figure S11A - Given that the oligomerization of N protein acts as a platform for viral RNA interactions 35 , we sought to examine the binding affinity of mutant and control N protein with viral RNA isolated from COVID-19 patient swabs. The RNA-binding activity of mutant and control N proteins was examined by pulled-down viral RNA through an in vitro RIP assay ( Figure 4B ), and our data revealed that the mutant N protein enriched significantly higher level of viral RNA compared to control protein ( Figure 4C ). This indicates a strong binding capability of mutant N proteins with viral RNA, which could potentially impact the essential roles of N protein at various stages of viral life cycle and its interaction with the host. According to the SIFT tool 29 , a substitution at position 204 from G to R in the N protein is predicted to affect functional properties ( Figure 4A ). Therefore, we decided to investigate how the two amino acids substitution (R203K and G204R) in the N protein impact its functional interaction with the host that could modulate viral pathogenesis and rewiring of host cell pathways and processes. HEK-293T cells (3 biological replicates, Supplementary Information) were used for affinity-purification followed by mass spectrometry analysis (AP-MS) to identify host proteins associated with the control and mutant N protein ( Figure S12 ). The majority (87%) of non-differentially interacting proteins overlapped with the previously reported 36 N protein interacting partners ( Figure S12D and Table S5 ). We identified 48 human proteins that displayed significant (adjusted p-value ≤ 0.05, and Log 2 fold change ≥ 1) differential interactions with the mutant and control N protein ( Figure 4D , Figure S12E , Table S6 ). Among these, 43 proteins . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint showed increased interaction and 5 proteins showed decreased interaction with the N mutant ( Figure 4D , Figure S12E ). Among the group with increased interaction, we identified many proteins associated with TOR and other signaling pathways (such as AKT1S1 and PIN1), proteins associated with the viral process, viral transcription, and negative regulation of RNA nuclear export (NUP153 and NUP98), and proteins involved in apoptotic and cell death processes (PAWR, ACIN1, and PDCD5) ( Figure 4D , Figure S12E ). We also identified proteins in the mutant condition that are linked with the immune system processes (PTMS), kinase activity (GCN1), and translation (e.g. MRPS36) ( Figure 4D , Figure S12E ). In the group with decreased interaction, we identified SNIP1 (NF-kappaB signaling), TMA16 (translation), and CSNK2B (casein kinase II) ( Figure 4D , Figure S12E ). Gene ontology analysis showed that the most enriched biological processes are associated with negative regulation of tRNA and ribosomal subunit export from the nucleus ( Figure 4E ). This finding suggests that the mutant virus may more efficiently inhibit and hijack the host translation to facilitate viral replication and pathogenesis. Further, many viruses can manipulate the host sumoylation process to enhance viral survival and pathogenesis 37 . By pathway enrichment analysis of differentially interacting proteins, we identified pathways associated with the sumoylation of host proteins and antiviral mechanisms ( Figure 4E ). In SARS-CoV, it has been shown that phosphorylation of the N protein is more prevalent during viral transcription and replication 38 is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint of N protein in the LKR region is critical for regulating both viral genome processing (transcription and replication) and nucleocapsid assembly 35, 40 . To further understand the functional relevance of KR mutation in the N protein, we performed phosphoproteomic analysis in control and mutant conditions. We consistently found that the serine 206 (S206) site, which is next to the KR mutation site ( Figure 4F) , is highly phosphorylated, specifically in the mutant N protein ( Figure 4G , Table S7 ). Notably, the phosphorylation level at serine 2 (S2) and other serine sites (S79, S176, and S180) within the LKR region did not change between mutant and control conditions ( Figure 4F ). Each node represents an enriched term and colored by its p-value. The size of each node corresponds to number of linked genes from the list. To understand whether the R203K/G204R mutations in the N gene affect host cell transcriptome, we transfected HEK293T cells with plasmids expressing the full-length N-control and N-mutant protein along with mock-transfection control. The transcriptome profile of Nmutant and N-control transfected cells displays a distinct pattern from the mock-control ( Figure S13A ). We identified 83 and 67 differentially expressed genes (DEGs) in the N-mutant and Ncontrol transfected cells, respectively, with adjusted p-value < 0.05 and log 2 fold-change ≥ 1 (Figure S13B-C and Table S8 ). Among the DEGs, numerous interferon, cytokine, and immunerelated genes are up-regulated, some of which are shown in Figure 5 (for complete list see Table S8 ). We found a robust overexpression of interferon-related genes in the N-mutant compared to N-control transfected cells ( Figure 5A -B) after adjusting for fold change ( Figure S13D ). Indeed, strong overexpression of interferon and chemokine related genes (Table S8) . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint Also, we found overexpression of other genes such as ACE2, STAT1 44 , and TMPRSS13 48 ( Figure 5A and Table S8 ) that are elevated in critical COVID-19 disease. Pathway enrichment analysis of the uniquely up-regulated genes in the N-mutant condition ( Figure 5C ) shows an overrepresentation of biological process pathways associated with response to the virus ( Figure 5D) . Similarly, all up-regulated genes were related to substantially enriched pathways, such as interferon-related response, cytokine production, and viral reproductive processes ( Figure S13E ). The enriched GO terms display an interconnected network highlighting the relationships between up-regulated overlapping genes sets in these pathways ( Figure 5D and Figure S13E ). Taken together, these results suggest that the R203K/G204R mutations in the N protein may enhance its function in provoking a hyper-expression of interferon-related genes that contribute to the cytokine storm in exacerbating COVID-19 pathogenesis. Figure 2C ). From estimates of viral genetic diversity and reproduction rate, we find that decreased diversity and reproduction rate coincides with imposed national curfews and is followed by an observed drop in reported COVID-19 cases ( Figure 1D ). . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint Our COVID-19 patient data allowed to us detect three SNPs -underlying the N protein R203K and G204R mutations -significantly associated with higher viral load. It is worth noting that two studies have found higher viral load has in infected patients to be associated with severity and mortality 49, 50 . Among our samples we initially observed an apparent association between the R203K/G204R mutations and mortality, however, the association was no longer statistically significant when correcting for sampling time. A dated phylogenetic approach suggests that the R203K/G204R mutations arose early in the SARS-CoV-2 pandemic -perhaps even as multiple independent events -and that the mutations entered Saudi Arabia during late January 2020, most likely through Italy. The N protein of SARS-CoV-2, an abundant viral protein within infected cells, serves multiple functions during viral infection, which besides RNA binding, oligomerization, and genome packaging, playing essential roles in viral transcription, replication, and translation 30, 51 . Also, the N protein can evade immune response and perturbs other host cellular processes such as translation, cell cycle, TGFβ signaling, and induction of apoptosis 52 to enhance virus survival. The critical functional regulatory hub within the N protein is a conserved serine-arginine (SR) rich-linker region (LKR), which is involved in RNA and protein binding 53 , oligomerization 33, 34 , and phospho-regulation 35, 40 . We show that the mutant N protein containing R203K and G204R changes has higher oligomerization and stronger viral RNA binding ability, suggesting a potential link of these mutations with efficient viral genome packaging. The R203K and G204R mutations are in close proximity to the recently reported RNA-mediated phase separation domain (aa 210-246) 42 that is involved in viral RNA packaging through phase separation. This domain was thought to enhance phase-separation also through protein-protein interactions 42 . Further studies are needed . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint to examine any definite link between KR mutation and phase-separation; however, the differential interaction of host proteins, as shown in our study could affect this process. Moreover, the functional activities of the N protein at different stages of the viral life cycle are regulated by phosphorylation-dependent physiochemical changes in the LKR region 40 As the COVID-19 pandemic is still ongoing, there is a need for novel therapeutic strategies to treat severe infections in patients. Our identified interaction pathways and inhibition of serine 206 phosphorylation could be used as potential targets for therapies. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint In conclusion, our results highlight the major influence of the R203K/G204R mutations on the essential properties and phosphorylation status of SARS-CoV-2 N protein that lead to increased host response and efficacy of viral infection. As part of the study, nasopharyngeal swab samples were collected in 1ml of TRIzol (Ambion, High quality SNPs were filtered using the filter expression: . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Consensus assembly sequences were deposited to GISAID (Table S1) 11 . To retrieve highconfidence SNPs assembled sequences were re-aligned against the Wuhan-Hu-1 reference sequence (NC_045512.2), and only positions in the sample sequences with unambiguous bases in a 7-nucleotide window centred around the SNP position were kept for further analysis. To generate the phylogeny of Saudi samples with a global context, a total of 308,012 global sequences were downloaded from GISAID on 31 December 2020, filtered and processed using Nextstrain pipeline 12 . Global sequences were grouped by country and sample collection month and 20 sequences per group were randomly sampled which resulted in 10,873 global representative sequences and 952 Saudi sequences. The phylogeny was constructed using IQ-TREE 64 , clades were assigned using Nextclade and internal node dates were inferred and sequences pruned using TreeTime 65 . Nextstrain protocol was followed for the above-mentioned steps. The resulting global phylogenetic tree was reduced to retain the branches that lead to Saudi leaf nodes and visualised using baltic library (https://github.com/evogytis/baltic). . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint Phylodynamic analyses use the same sequence subset used in the full phylogenetic analysis, extracted from the GISAID SARSCoV-2 database 11 . Wrapper functions for the importation date estimates and skygrowth model are provided in the sarscov2 R package as 'compute_timports' and 'skygrowth1' respectively (github.com/emvolz-phylodynamics/sarscov2Rutils). Sequences corresponding to each Nextstrain 12 clade were extracted using the Nextstrain_clade parameter in the GISAID metadata table. A subset of 500 international sequences were select for each clade based on Tamura Nei 93 distance with tn93 (github.com/veg/tn93) and stratified over time 66 To identify import events that resulted in new introductions into Saudi Arabia, 25,198 sequences were subsampled from 590K global sequences available on GISAID on February 24th 2021. Samples with closer genetic distance to Saudi Samples were preferred. The phylogeny was constructed using IQ-TREE 64 , internal nodes dates and possible country for internal nodes were inferred using TreeTime 65 . Nextstrain protocol was followed for the above-mentioned steps 12 . In . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. We modelled growth rate and effective population size over time on these phylogenies using the R package skygrowth 16 . Skygrowth is a non-parametric Bayesian approach which applies a stochastic process on estimates of growth rate and effective population size. The model included mean-centred, unit variance estimates of travel rates from google mobility data (google.com/covid19/mobility/) as a covariate (transit stations percent change from baseline), 60 timesteps and a tau (precision) value corresponding to a 1% change in growth per week. The growth rate output was converted to an estimate of R over time using an infectious period of 9.5 days 70 . Origin of R203K/G204R SNPs A total of 590K samples submitted to GISAID until February 24 were downloaded and SNPs indentifed by mapping against the Wuhan reference using minimap2 71 . The variants were queried to count the distribution of triplets among various Nextstrain clades ( Figure S7 ). To identify if there are lineages of triplet SNPs in clades other than 20B, a phylogenetic tree was . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint constructed by including all R203K/G204R samples found in other clades outside 20B and its subclades ( Figure S4 ). As it was already evident that 20B and its subclades contains lineages of R203K/G204R samples, subsamples from 20B and its subclades were sufficient to obtain a total of 16,386 samples. Statistical analyses were performed with the statistical software R version 4.0.3 72 Table S9 . Cell culture and transfection HEK293T (ATCC; CRL-3216) cells were grown in Dulbecco's modified Eagle's medium (DMEM) (4.5 g/l d-glucose and Glutamax, 1 mM sodium pyruvate) (GIBCO) and 10% fetal bovine serum (FBS; GIBCO) with penicillin-streptomycin supplement, according to standard protocols (culture condition 37 °C and 5% CO2). Transfection of ten million cells per 15-cm dish . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. protease and phosphatase inhibitor cocktails) and then proceed with on-bead digestion. The onbead digestion was carried out as described before 36 . For affinity confirmation, bound proteins were eluted using buffer BXT (IBA Lifesciences) and after running on SDS-PAGE were subjected to silver staining and western-blot using anti-strep-II antibody (ab76949). To purify clean 2xStrep-tagged N protein (mutant and control), we applied stringent washing and double elution strategy. The MS analysis was performed as described previously 73, 74 with slight modifications. For mass spectrometry analysis an Orbitrap Fusion mass spectrometer (MS) (Lumos, Thermo Fisher . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint Scientific) was used in data-dependent acquisition (DDA) mode. For injection, 0.5 µg peptide mixture was used and desalting was performed for 5 minutes in 0.1% FA in water. The gradient and all other steps were essentially the same as described 73 . Protein identification analysis from the raw mass spectrometry data was performed using the Maxquant software (version 1.5.3.30) 75 as described 73 . For phosphorylated peptides, we used Maxquant label-free quantification (LFQ) 75 . The analysis and quantification of phosphorylated peptides was performed according to published protocol 76 . The normalized LFQ data were processed for statistical analysis on the LFQ-Analyst a webbased tool 77 to performed pair-wise comparison between mutant and control N protein AP-MS data. The significant differentially changed proteins between mutant and control conditions were identified. The threshold cut-off of adjusted p-value <= 0.05, and Log fold change >= 1 were used. Among the replicates, outliers were removed based on correlation and PCA analysis. The GO enrichment analysis was performed on the LFQ-Analyst 77 . Bis(sulfosuccinimidyl) suberate (BS3, Thermo Scientific Pierce) was used for cross-linking of control and mutant N protein to analyse the oligomerization properties. The experiment was performed as reported previously 78 . RNA-sequencing and differential gene expression analysis HEK293T cells were transfected with plasmids expressing the full-length N-control and Nmutant protein along with mock control. After 48-hour cells were harvested in Trizol and total . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint RNA was isolated using Zymo-RNA Direct-Zol kit (Zymo, USA) according to the manufacture's instruction. The concentration of RNA was measured by Qubit (Invitrogen), and RNA integrity was determined by Bioanalyzer 2100 system (Agilent Technologies, CA, USA). The RNA was then subjected to library preparation using Ribozero-plus kit (Illumina). The libraries were sequenced on NovaSeq 6000 platform (Illumina, USA) with 150 bp paired-end reads. The raw reads from HEK293T RNA-sequencing were processed and trimmed using trimmomatic 60 and mapped to annotated ENSEMBL transcripts from the human genome (hg19) 79,80 using kallisto 81 . Differential expression analysis was performed after normalization using EdgeR integrated in the NetworkAnalyst 82 . GO biological process and pathway enrichment analyses on up-regulated genes were performed using NetworkAnalyst 82 . We are sincerely grateful to all hospital members for providing samples and collating metadata in such an unprecedented pandemic, along with the MOH and ethical committee, which rendered it permissible. We thank the KAUST Rapid Research Response Team (R3T) under the Vice President -Research (VPR) office in KAUST for generous financial support. We also thank Erik Talley from KAUST Health Safety and Environment (HSE) and Hani Bukhari from KAUST Security for providing timely logistical support for samples transport during COVID-19 Curfew restrictions in the Kingdom. We extend our thanks and appreciation to GDRS director, PH. Athari Alotaibi (General Director for Research and Studies, MOH) for her vigorous facilitation of the research project, and Mohammad Fawzi (General Directorate of Health Affairs) for his help with the metadata collection. We also deeply thank Dr. Wael Hamzah Motair, Dr. Nader Hamzah Motair, Dr. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. We gratefully acknowledge all of the authors from the originating laboratories responsible for obtaining the specimens and the submitting laboratories where genetic sequence data were generated and publicly shared via the GISAID Initiative, on which was partially used for additional support for some of the conclusions drawn in this study. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Supplementary Information is available for this paper. Assembled virus genomes are available at GISAID (Table S1 ). Reads from RNA-Seq analysis of transfected HEK293T cells have been uploaded to European Nucleotide Archive (https://www.ebi.ac.uk/ena/) under the Study accession number PRJEB44716. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 10, 2021. ; https://doi.org/10.1101/2021.05.06.21256706 doi: medRxiv preprint Coronavirus disease (COVID-19) Weekly Epidemiological Update and Weekly Operational Update, An interactive web-based dashboard to track COVID-19 in real time COVID-19 Dashboard COVID-19: preparing for superspreader potential among Umrah pilgrims to Saudi Arabia Tale of three seeding patterns of SARS-CoV-2 in Saudi Arabia Estimation of Coronavirus Disease 2019 (COVID-19) Burden and Potential for International Dissemination of Infection From Iran Saudi Arabia announces first case of coronavirus Genomic determinants of pathogenicity in SARS-CoV-2 and other human coronaviruses Genomic Epidemiology of SARS-CoV-2 in Guangdong Province Data, disease and diplomacy: GISAID's innovative contribution to global health Global initiative on sharing all influenza data -from vision to reality Nextstrain: real-time tracking of pathogen evolution Evaluating the Effects of SARS-CoV-2 Spike Mutation D614G on Transmissibility and Pathogenicity Increased mortality in community-tested cases of SARS-CoV-2 lineage B.1.1.7. Nature Genomic monitoring of SARS-CoV-2 uncovers an Nsp1 deletion variant that modulates type I interferon response Modeling the Growth and Decline of Pathogen Effective Population Size Provides Insight into Epidemic Dynamics and Drivers of Antimicrobial Resistance Three adjacent nucleotide changes spanning two residues in SARS-CoV-2 nucleoprotein: possible homologous recombination from the transcription-regulating sequence. bioRxiv Phylogenetic analysis of SARS-CoV-2 data is difficult Stability of SARS-CoV-2 phylogenies An integrated national scale SARS-CoV-2 genomic surveillance network Novel Coronavirus Is Undergoing Active Recombination No detectable signal for ongoing genetic recombination in SARS-CoV-2. bioRxiv Effects of SARS-CoV-2 mutations on protein structures and intraviral protein-protein interactions A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology Recombinant SARS-CoV-2 genomes involving lineage B.1.1.7 in the UK Generalized Additive Models: An Introduction with R, 2 edition Tracking Changes in SARS-CoV-2 Spike: Evidence that D614G Increases Infectivity of the COVID-19 Virus SIFT: predicting amino acid changes that affect protein function The coronavirus nucleocapsid is a multifunctional protein Evolutionary dynamics of SARS-CoV-2 nucleocapsid protein and its consequences A genetic barcode of SARS-CoV-2 for monitoring global distribution of different clades during the COVID-19 pandemic Analysis of multimerization of the SARS coronavirus nucleocapsid protein Transient Oligomerization of the SARS-CoV N Protein -Implication for Virus Ribonucleoprotein Packaging Characterization of SARS-CoV-2 N protein reveals multiple functional consequences of the C-terminal domain A SARS-CoV-2 protein interaction map reveals targets for drug repurposing Viral manipulation of the cellular sumoylation machinery Nucleocapsid Phosphorylation and RNA Helicase DDX1 Recruitment Enables Coronavirus Transition from Discontinuous to Continuous Transcription Glycogen Synthase Kinase-3 Regulates the Phosphorylation of Severe Acute Respiratory Syndrome Coronavirus Nucleocapsid Protein and Viral Replication Phosphoregulation of Phase Separation by the SARS-CoV-2 N Protein Suggests a Biophysical Basis for its Dual Functions Nucleocapsid protein of SARS-CoV-2 phase separates into RNA-rich polymerase-containing condensates The SARS-CoV-2 nucleocapsid phosphoprotein forms mutually exclusive condensates with RNA and the membrane-associated M protein Transcriptional profiling of leukocytes in critically ill COVID19 patients: implications for interferon response and coagulation Host transcriptomic profiling of COVID-19 patients with mild, moderate, and severe clinical outcomes Two distinct immunopathological profiles in autopsy lungs of COVID-19 In vivo antiviral host transcriptional response to SARS-CoV-2 by viral load, sex, and age Severity of SARS-CoV-2 infection as a function of the interferon landscape across the respiratory tract of COVID-19 patients. bioRxiv TMPRSS11D and TMPRSS13 Activate the SARS-CoV-2 Spike Protein SARS-CoV-2 viral load is associated with increased disease severity and mortality SARS-CoV-2 viral load predicts COVID-19 mortality The SARS coronavirus nucleocapsid protein--forms and functions Molecular Biology of the SARS-Coronavirus View from an mRNP: The Roles of SR Proteins in Assembly, Maturation and Turnover The Global Phosphorylation Landscape of SARS-CoV-2 Infection The Multifarious Role of 14-3-3 Family of Proteins in Viral Replication The Coronavirus Nucleocapsid Protein Is Dynamically Associated with the Replication-Transcription Complexes Mass spectroscopic characterization of the coronavirus infectious bronchitis virus nucleoprotein and elucidation of the role of phosphorylation in RNA binding by using surface plasmon resonance Phosphorylation of the arginine/serine dipeptiderich motif of the severe acute respiratory syndrome coronavirus nucleocapsid protein modulates its multimerization, translation inhibitory activity and cellular localization Determination of host proteins composing the microenvironment of coronavirus replicase complexes by proximity-labeling Trimmomatic: a flexible trimmer for Illumina sequence data Fast and accurate long-read alignment with Burrows-Wheeler transform The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies Maximum-likelihood phylodynamic analysis Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees Dating of the human-ape splitting by a molecular clock of mitochondrial DNA Scalable relaxed clock phylogenetic dating phylogenetic analysis in R Clinical characteristics of 24 asymptomatic infections with COVID-19 screened among close contacts in Nanjing Minimap2: pairwise alignment for nucleotide sequences R: A language and environment for statistical computing. R Foundation for Statistical Computing Arabidopsis proteome and the mass spectral assay library Ubiquitin ligases HUWE1 and NEDD4 cooperatively control signal-dependent PRC2-Ezh1 alpha/beta-mediated adaptive stress response pathway in skeletal muscle cells MaxQuant enables high peptide identification rates, individualized p.p.b.-range mass accuracies and proteome-wide protein quantification Glucose-regulated phosphorylation of TET2 by AMPK reveals a pathway linking diabetes to cancer LFQ-Analyst: An Easy-To-Use Interactive Web Platform To Analyze and Visualize Label-Free Proteomics Data Preprocessed with MaxQuant Biochemical characterization of SARS-CoV-2 nucleocapsid protein The Ensembl gene annotation system Near-optimal probabilistic RNA-seq quantification NetworkAnalyst 3.0: a visual analytics platform for comprehensive gene expression profiling and meta-analysis