SARS-CoV-2 genomes from Saudi Arabia implicate nucleocapsid mutations in host response and increased viral load

Monitoring SARS-CoV-2 spread and evolution through genome sequencing is essential in handling the COVID-19 pandemic. Here, we sequenced 892 SARS-CoV-2 genomes collected from patients in Saudi Arabia from March to August 2020. We show that two consecutive mutations (R203K/G204R) in the 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 increased interaction of GSK3A kinase simultaneously with hyper-phosphorylation of the adjacent serine site (S206) in the mutant N protein. Furthermore, the host cell transcriptome analysis suggests that the mutant N protein produces dysregulated interferon response genes. Here, we provide crucial information in linking the R203K/G204R mutations in the N protein to modulations of host-virus interactions and underline the potential of the nucleocapsid protein as a drug target during infection.

T he emergence of novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes the respiratory coronavirus disease 2019 , 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 6457 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 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 2nd, 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 . As of October 2021, over four and a half million complete and high coverage genomes are accessible on GISAID 10,11 , which aids immensely in tracking the 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 (Fig. 1a, b, Supplementary Data 1, Supplementary Table S1, and  Supplementary Table S2). Samples were mainly collected early in the pandemic and during the first wave of reported cases in Saudi Arabia (Fig. 1c). 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 the 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 a specific set of mutation events in the N gene to viral loads in a diverse population of COVID-19 patients in Saudi Arabia (Supplementary Fig. S1). Finally, we experimentally show the functional impact of these mutations in the N protein on the virus' interactions with the host.

Results
SNP calling and phylodynamics of SARS-CoV-2 samples from Saudi Arabia. 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 (Supplementary Table S1). Data on comorbidities were available for 689 patients with diabetes (39%) and hypertension (35%) being the most abundant (Supplementary Table S2). Patient outcome data were available for 850 samples, and 199 patients (23%) died during hospitalization (Supplementary Table S1).
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 SARS-CoV-2 Wuhan-Hu-1 isolate reference (GenBank accession: NC_045512) (Supplementary Fig. S2). The observed numbers of SNPs relative to the reference sequence are in general lower than the numbers observed in global samples, but with the exception of a period from mid-June to late July, the average number of SNPs in Saudi samples is within one standard deviation of samples deposited in GISAID (Supplementary Fig. S3). We further detected 41 indels of which 26 reside in coding regions (Supplementary Table S3). 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 (Supplementary Fig. S2). These include the Spike protein D614G (A23403G) and three consecutive SNPs (G28881A, G28882A, and G28883C) causing the R203K and G204R changes in the nucleocapsid protein. Together with all sequences from Saudi Arabia available on GISAID on December 31st 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 (Fig. 1d) shows a downward trend in the effective reproduction number (R(t)) over time with the timely introduction and maintenance of effective non-pharmaceutical interventions by the Saudi Ministry of Health. Despite the bounds of R(t) estimate including one for much of the study period, we can infer an estimated decrease in R(t) over time until the lifting of restrictions in late June (Fig. 1d). By investigating the 1.6 million individual traces from the Skygrowth analysis we can infer a date of 27th April 2020 as the first point where over 50% of collected traces resulted in an R(t) estimate below one, suggesting an epidemic in decline for the first time.
The effective population size (Ne) represents the relative diversity of the sequences collected in Saudi Arabia over the course of the outbreak (Fig. 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 five major Nextstrain clades 12 , 19A-B and 20A-C (Fig. 2a). The approximate relationships between Nextstrain clades and Pangolin lineages 17 are listed in Supplementary Table S4). The phylogeny highlighted the clade 20A that all carried the nucleocapsid (N) protein R203K/G204R mutations 18 with high incidences of ICU hospitalizations. These samples were predominantly coming from Jeddah. Through timescaled phylogenies dates of importation events were then estimated for each clade. The estimates of importation date with phangorn (see "Methods") are early, with clade 19B estimated to have been imported with the highest probability in late February 2020 and the other reported clades most likely imported in March (Fig. 2b). As this method cannot directly estimate the timing of an importation event, importations were estimated to occur at the midpoint of a branch along which a state change occurs between the internal and external samples. As estimated branches are short due to the inclusion of closely related external sequences this should not significantly alter the estimated importation dates but does make an assumption about the timing of the movement event which could be overcome with more complex and computationally intensive phylogeographic methods and additional population movement data.
Origin of R203K/G204R SNPs. 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 ( Supplementary Fig. 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 constructing a reliable and robust phylogeny is problematic 19 , 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 20,21 . Additionally, the placement of individual virus genomes may be hampered by systematic errors, homoplasies, potential recombination, or coinfection of multiple virus strains 19,20,[22][23][24][25][26][27] . 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 ( Supplementary Fig. S4). Within our sampling window we observe an apparent transient increase in the frequency of R203K/G204R SNPs (Fig. 3a) in accordance with earlier observations 18,28 . In the global data, the peak in R203K/ G204R frequency is slightly delayed compared to samples from Saudi Arabia. The initial global peak is observed in July 2020 followed by a decline until the fall of 2020, where the R203K/ G204R SNPs once again increased along with the Spike protein Y501N mutation in the B1.1.17 lineage 17 (Fig. 3a).
A mutant form of the nucleocapsid (N) protein associated with higher viral loads in COVID-19 patients in Saudi Arabia. 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 (Fig. 3b, Supplementary Fig. S5). 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 ( Supplementary  Fig. S6). In our samples, no other SNPs co-occur with the R203K/ G204R SNPs ( Supplementary Fig. S7). 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 (Supplementary Table S1).
Using multivariable regression, we next evaluated the effect of the R203K/G204R SNPs on mortality, severity, and viral load in our COVID-19 patient samples for which a 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 logistic linear model using R203K/G204R SNPs as a covariate and adjusting by sex, age, comorbidities, hospital, and other SNPs. We included 12 additional SNPs (C241T, C1191T, C3037T, G10427A, C14408T, C15352T, C18877T, A23403G, G25563T, C26735T, T27484C, and C28139T) that co-occurred with the R203K/G204R SNPs in at least five samples in the model. The A23403G mutation results in the Spike protein D614G SNP that is associated with higher viral load 13 . Then, to adjust for confounding effects from time, we also fitted models that included time. In the models, age and time were included using smoothing splines to allow for potential nonlinear relationships 29 .  Using first a logistic regression model that did not include time, 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.18, 95% CI 0.22-2.13 (Supplementary Table S5). A positive significant association was also observed for the C14408T SNP, and a negative association for the C241T SNP (Supplementary Table S5).
In the time-adjusted model, the log-odds for the R203K/G204R SNPs increased to 1.38, 95% CI 0.28-2.48 (Supplementary  Table S6). In this model, the C241T SNP again displayed a significant negative association, and a positive association was now observed for the C1887T SNP (Supplementary Table S6).
The relationship between mortality and R203K/G204R SNPs was positive and statistically significant in the model that did not include time with log-odds equal to 1.04, 95% 0.16-1.92. No significant association was observed for other SNPs (Supplementary Table S7). However, after adjusting for time as a variable, there was no longer any association between R203K/G204R SNPs and mortality (log-odds: 0.58, 95% CI −0.41-1.56) (Supplementary Table S8). 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 ( Supplementary Fig. S8). One could speculate if temporal shifts in clinical settings in certain cases have impacted the balance between the severity of infection and mortality. Unfortunately, our dataset 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. As two different kits were used for the qPCR reactions (see Methods), we fitted adjusted models that besides sex, age, comorbidities, hospital, and time, included qPCR kits and the above-mentioned SNPs as covariates. From this adjusted 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.33 units (95% CI 0.72-1.93) ( Supplementary Table S9). Similarly, the model showed a positive significant association between the SNPs A23403G (Spike protein D614G) and C26735T SNPs and log 10 (viral copy number), the former being consistent with earlier reports 13,30 . A significant negative association was found for the C3037T C14408T, and G25563T SNPs (Supplementary Table S9). The positive and statistically significant association of R203K/G204R SNPs with higher viral load in critical COVID-19 patients (Fig. 3c) hence suggests their functional implications during viral infection.
N mutant protein has high oligomerization potential and RNA-binding affinity. The SARS-CoV-2 N protein binds the viral RNA genome and is central to viral replication 31 . Protein structure predictions have shown that the R203K/G204R mutations result in significant changes in protein structure 28 , theoretically destabilizing the N structure 32 , and potentially enhancing the protein's ability to bind RNA and alter its response to serine  33 . 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) (Fig. 4a), known to be involved in the oligomerization of N proteins 34,35 . Protein cross-linking shows that N mutant protein (with the R203K/G204R mutations) has higher oligomerization potential compared to the control N protein (without the changed amino acids) at low protein concentration ( Supplementary  Fig. S9).
Given that the oligomerization of N protein acts as a platform for viral RNA interactions 36 , we sought to examine the binding affinity of mutant and control N protein with viral RNA isolated ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28287-8 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 (Fig. 4b), and our data revealed that the mutant N protein enriched significantly higher level of viral RNA compared to control protein (Fig. 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.
The R203K/G204R mutations in the N protein affect its interaction with host proteins. According to the SIFT tool 37 , a substitution at position 204 from G to R in the N protein is predicted to affect functional properties (Fig. 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 (three biological replicates) were used for affinity purification followed by mass spectrometry analysis (AP-MS) to identify host proteins associated with the control and mutant N protein ( Supplementary Fig. S10). Protein identification was performed using MaxQuant software 38 . The compilation of the identified protein groups in mock and N protein (N control and mutant) AP-MS is presented in Supplementary Data 2 (Raw data are available via ProteomeXchange with identifier PXD027168). The majority (62%) of previously reported 39 N protein interacting partners overlapped with the identified unadjusted proteins list ( Supplementary Fig. S10). We identified 43 human proteins that displayed significant (adjusted p-value ≤ 0.05, and Log 2 foldchange ≥ 1) differential interactions with the mutant and control N protein (Fig. 4d, Supplementary Fig. S10, and see Supplementary Data 3 for both significant and non-significant differentially interacting protein list). Among these, 42 proteins showed increased interaction and one protein (PRPF19) showed decreased interaction with the N mutant ( Fig. 4d and Supplementary Fig. S10). 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 and ACIN1) ( Fig. 4d and Supplementary Fig. S10). We also identified proteins in the mutant condition that are linked with the immune processes and translation ( Fig. 4d and Supplementary Fig. S10). Gene ontology analysis showed that the most enriched biological processes are associated with negative regulation of tRNA and ribosomal subunit export from the nucleus (Fig. 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 40 . By pathway enrichment analysis of differentially interacting proteins, we identified pathways associated with sumoylation and antiviral mechanisms (Fig. 4e).
Serine 206 (S206) displays hyper-phosphorylation in the mutant N protein. In SARS-CoV, it has been shown that phosphorylation of the N protein is more prevalent during viral transcription and replication 41 and inhibition of phosphorylation diminishes viral titer and cytopathogenic effects 42 . Recent elegant studies elaborated the role of N protein phosphorylation in modulating RNA binding and phase separation in SARS-CoV-2 36,[43][44][45] . Thus, phosphorylation of N protein in the LKR region is critical for regulating both viral genome processing (transcription and replication) and nucleocapsid assembly 36,43 . 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 (Fig. 4f), is highly phosphorylated, specifically in the mutant N protein ( Fig. 4g and Supplementary Data 4). Notably, the phosphorylation level at serine 2 (S2) and other serine sites (S79, S176, and S180) within the LKR region is detectable in both mutant and control conditions (Supplementary Data 4 shows only N protein phospho(STY) sites data. Complete raw data are available with the identifier PXD027168).
The N mutant (R203K/G204R) induces overexpression of interferon-related genes in transfected host cells. To understand whether the R203K/G204R mutations in the N gene affect host cell transcriptome, we transfected Calu-3 cells (4 biological replicates) with plasmids expressing the full-length N-control and N-mutant protein along with mock-transfection control. The transcriptome profile of N-mutant and N-control transfected cells displays a distinct pattern from the mock control ( Supplementary  Fig. S11). We identified 144 and 153 differentially expressed (DE) genes in the N-control and N-mutant transfected cells, respectively, with adjusted p-value < 0.05 and log 2 fold-change ≥1 ( Supplementary Fig. S11 and Supplementary Data 5). Among the DE genes, numerous interferon, cytokine, and immune-related genes are up-regulated, some of which are shown in Fig. 5 (for complete list see Supplementary Data 5). We found a robust overexpression of interferon-related genes in the N-mutant compared to N-control transfected cells (Fig. 5a-b) after  Fig. S11). Indeed, strong overexpression of interferon and chemokine-related genes (Supplementary Data 5) were reported in critical COVID-19 patients 46,47 . Recent reports further indicate a link between increased expression of interferon-related genes and higher viral load in severe COVID-19 patients [48][49][50] . Also, we found overexpression of other genes such as ACE2, STAT1 47 , and TMPRSS13 51 (Fig. 5a and Supplementary Data 5) that are elevated in critical COVID-19 disease.
Pathway enrichment analysis (top 15 pathways based on p-value and FDR) of the up-regulated genes (Fig. 5c) shows an overrepresentation of biological processes associated with response to the virus ( Fig. 5d and Supplementary Data 6). Similarly, all DE genes were related to substantially enriched pathways, such as interferonrelated response, cytokine production, and viral reproductive processes ( Supplementary Fig. S11). The enriched GO terms display an interconnected network highlighting the relationships between up-regulated overlapping genes sets in these pathways ( Fig. 5d and Supplementary Fig. S11). 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. influx of SARS-CoV-2 lineages into the Kingdom. The earliest estimated importation dates point to an entry during the early stages of the pandemic (Fig. 2b). 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 (Fig. 1c, d).

Discussion
Our COVID-19 patient data allowed to us detect three SNPsunderlying 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 52,53 . 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. The association of N protein R203K and G204R mutations with higher viral load persists after adjusting for time but not with mortality, suggesting that the mortality rate of severe infections may be influenced by other factors such as changes in treatment regimes as well as the complexities of host response to SARS-CoV-2 infection. Unfortunately, the available data do not allow a further assessment of this.
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 31,54 . 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 55 to enhance virus survival. The critical functional regulatory hub within the N protein is a conserved serinearginine (SR) rich-linker region (LKR), which is involved in RNA and protein binding 56 , oligomerization 34,35 , and phosphoregulation 36,43 .
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) 45 that is involved in viral RNA packaging through phase separation. This domain was thought to enhance phase separation also through protein-protein interactions 45 . Further studies are needed 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 phosphorylationdependent physiochemical changes in the LKR region 43 . Although all individual phosphorylation sites may not be functionally important 33,57 , the specific enhancement of phosphorylation at serine 206 in the mutant N protein shown in this study hints at its functional significance. The serine 206 can form a phosphorylation-dependent binding site for protein 14-3-3, involved in cell cycle regulatory pathways regulating human and virus protein expression 58 . Multiple lines of evidence show that N protein phosphorylation is critical for its dynamic localization and function at replication-transcription complexes (RTC), where it promotes viral RNA transcription and translation by recruiting cellular factors [41][42][43][59][60][61][62] . The enrichment of glycogen synthase kinase 3 A (GSK3A) with the mutant N protein, could specifically phosphorylate serine 206 in the R203K/G204R mutation background. GSK3 was shown to be a key regulator of SARS-CoV replication due to its ability to phosphorylate N protein 42 . Phosphorylation of serine 206 acts as priming site for initiating a cascade of GSK-3 phosphorylation events 42,43 . Also, GSK3 inhibition dramatically reduces the production of viral particles and the cytopathic effect in SARS-CoV-infected cells 42 . The detection of hyperphosphorylation at serine 206 only in the N-mutant could also point toward the phosphorylation timing difference or more phosphorylation stability at this site in the mutant background. Finally, our analysis of the transcriptome in transfected cells suggests that the mutant N protein induce overexpression of interferon-related genes that can aggravate the viral infection by inducing cytokine storm.
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.
In conclusion, our results suggest a major influence of the R203K/G204R mutations on the essential properties and phosphorylation status of SARS-CoV-2 N protein that may lead to increased host response and heightened efficacy of viral infection.

Methods
Sample collection. As part of the study, nasopharyngeal swab samples were collected in 1 ml of TRIzol (Ambion, USA) from 892 COVID-19 patients with various grades of clinical disease manifestations-consisting of severe, mild and asymptomatic symptoms. The anonymized samples were amassed from 8 hospitals and one quarantine hotel located in Madinah, Makkah, Jeddah, and Riyadh (Supplementary Table S1). Patient metadata in the form of age, sex, comorbidities, ICU submission, and mortality were provided by the hospitals (Supplementary  Table S2) and used for statistical analysis. Ethical approvals were obtained from the Institutional review board of the Ministry of Health in the Makkah region with the numbers H-02-K-076-0420-285 and H-02-K-076-0320-279, as well as the Institutional review board of Dr. Sulaiman Al Habib Hospital number RC20.06.88 for samples from Riyadh and the Eastern regions respectively. RNA Isolation. RNA was extracted using the Direct-Zol RNA Miniprep kit (Zymo Research, USA) following the manufacturer's instructions, along with several optimization steps to improve the quality and quantity of RNA from clinical samples. The optimization included extending the TRIzol incubation period, and the addition of chloroform during initial lysis step to obtain the aqueous RNA layer. The quality control of purified RNA was performed using Broad Range Qubit kit (Thermo Fisher, USA) and RNA 6000 Nano LabChip kit (Agilent, USA) respectively. RT-PCR was conducted using the one-step Super Script III with Platinum Taq DNA Polymerase (Thermo Fisher, USA) and TaqPath COVID-19 kit (Applied Biosystems, USA) on the QuantStudio 3 Real-Time PCR instrument (Applied Biosystems, USA) and 7900 HT ABI machine. The primers and probes used were targeting two regions in the nucleocapsid gene (N1 and N2) in the viral genome following the Centre for Disease Control and prevention diagnostic panel, along with primers and probe for human RNase P gene (CDC; fda.gov/media/ 134922/download) (Supplementary Table S10). Samples were considered COVID positive once the cycle threshold (Ct) values for both N1 and N2 regions were less than 40. For amplicon seq purposes, the samples chosen were of Ct less than 35 to ensure successful genome assembly in order to upload on GISAID.
Consensus assembly sequences were deposited to GISAID (Supplementary Data 1) 11 . To retrieve high-confidence SNPs assembled sequences were re-aligned Phylogenetic 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 (v2.0.5) 67 , clades were assigned using Nextclade and internal node dates were inferred and sequences pruned using TreeTime (v0.7.5) 68 . 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 visualized using baltic library (https://github.com/evogytis/baltic).
Phylodynamic analysis. 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 'sky-growth1' respectively (https://github.com/JorgensenD/sarscov2Rutils) 69 .
Importation date estimates for Nextstrain clades. Importation rate estimates were carried out using all available sequence data for Saudi Arabia deposited on GISAID 11 up to 31 December 2020, including the sequences described in this paper. Sequences were grouped by Nextstrain clade for analysis using the Next-strain_clade parameter 12 in the GISAID metadata table. Additional international sequences were selected for each of the included Nextstrain clade based on Tamura Nei 93 distance with the C program tn93 (v1.0.6) (github.com/veg/tn93) 70 . Five hundred sequences were selected from available closely related sequences in a time stratified manner, taking every N/500th sequence from the set of N sequences arranged by date, rounded to the nearest integer.
For each Nextstrain clade a maximum-likelihood phylogeny was produced with IQTree (v1.6.12) with an HKY substitution model 67,71 . These trees were dated using the R package treedater (v0.5.0) after collapsing short branch lengths and resolving polytomies randomly fifteen times for each clade with the functions di2multi and multi2di from the ape R package (v5.5) 72,73 . A strict molecular clock was used when estimating dated phylogenies, constrained between 0.0009 and 0.0015 substitutions per site per year 74 . The state of each internal node in the phylogeny was reconstructed by maximum parsimony with the R package phangorn (v2.7.0) 75 . As this method cannot directly estimate the timing of an importation event, importations were estimated to occur at the midpoint of a branch along which a state change occurs between the internal and external samples.
The probability density of importation events over time into Saudi Arabia by cluster is presented in Fig. 2b where a gaussian kernel density estimator is used with a bandwidth of 5.51 chosen using the Silverman rule of thumb, implemented in the geom_density function of the ggplot2 R package (v3.3.5) 76,77 .
Skygrowth model. The assembled sequences, together with all other reported sequences for Saudi Arabia available on GISAID on 31 December 2020, were used to construct estimates of the effective population size and growth rate of SARS-CoV2 in Saudi Arabia over the course of the first wave of the epidemic (March to September 2020) 11 . The collected sequences were used to produce a maximum likelihood phylogeny with an HKY substitution model in IQTree (v1.6.12) 67,71 . A set of 1000 bootstrap pseudo-replicate trees were produced with the ultrafast bootstrap approximation 78 . For each bootstrap phylogeny, branches with length less than 10 −5 were collapsed and polytomies were resolved randomly using the di2multi and multi2di functions in the ape R package (v5.5) 72 . These dichotomous trees were produced for subsequent molecular clock analysis and coalescent analysis. The set of 1000 initial bootstrap trees and 1000 additional phylogenies with randomly resolved short branches were dated using the reported collection dates for the sequenced samples as the date of the corresponding tip in each phylogeny. A strict molecular clock was used when estimating dated phylogenies, constrained between 0.0009 and 0.0015 substitutions per site per year with the R package treedater (v0.5.0) 73,74 .
The R package Skygrowth (v0.3.1) was used with these phylogenies to estimate the growth rate and effective population size of SARS-CoV2 in Saudi Arabia over time 16 . Skygrowth is a Bayesian non-parametric model of effective population size, with the primary difference to other skyline smoothing methods being that the first-order stochastic process is defined in terms of the growth rate of the effective population size (Ne) rather than Ne itself. The growth rate of Ne often has a simple relationship with the growth rate of the epidemic from which samples are collected.
Here, under the assumption of a simple SIR transmission process, we return the effective reproduction number as in Eq. (1) where _ Ne ðtÞ is the derivative of Ne ðtÞ , ð _ Ne ðtÞ =Ne ðtÞ Þ is the growth rate of Ne, and ψ is the mean generation time.
The model included 35 timesteps and an exponential prior on the smoothing parameter tau (precision) corresponding to a 1% change in growth per week. The growth rate output was converted to an estimate of R over time using an assumed generation time ψ of 9.5 days 79 . The model is fitted from March 5th to ensure exponential growth is established in the population prior to model fitting as we cannot estimate the growth rate or effective reproduction number when virus circulation is driven purely by importation events.
Although phylogenetic methods are sensitive to sampling rate changes over time, the coalescent method used in Skygrowth is relatively robust to heterogenous sampling through time, but can still be biased by unequal sampling in space or risk groups. Computation of reproduction numbers is premised on equivalence of the growth rate of the epidemic and the growth rate of the effective population size, which does not hold when the transmission rate is highly variable.
Origin of R203K/G204R SNPs. A total of 590 K samples submitted to GISAID until February 24 were downloaded and SNPs indentifed by mapping against the Wuhan-Hu-1 reference sequence (NC_045512.2) using minimap2 (v2.17) 80 . The variants were queried to count the distribution of triplets among various Nextstrain clades ( Supplementary Fig. S6). To identify if there are lineages of triplet SNPs in clades other than 20B, a phylogenetic tree was constructed by including all R203K/ G204R samples found in other clades outside 20B and its subclades (Supplementary Fig. 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.
Genome-wide associations. Tests were carried out as a conventional case-control setup, and for each SNP a contingency table was constructed as (deceased patients with SNP/ deceased patients without SNP)/(non-deceased patients with SNP/ nondeceased patients without SNP). P-values were calculated using Fisher's exact twotailed test.
Allele frequency estimate. Duplicate reads were removed from the mapped short Illumina sequence reads using picard tools' MarkDuplicates function (v2.20.4, Broad Institute, GitHub repository. http://broadinstitute.github.io/picard/). Frequencies were collected directly from the output of samtools (v1.8) mpileup 66 . Distal read positions were excluded and only read mapping and base calling qualities of at least 30 were considered.
Multivariable regression analysis. Statistical analyses were performed with the statistical software R version 4.0.3 81 and the R package mgcv (v1.8.33). Models for mortality and severity used 892 observations (all data).
Affinity purification and on-bead digestion. Cell lysis and affinity purification with MagStrep beads (IBA Lifesciences) were manually performed according to the published protocol 39 with minor modifications. Briefly, after transfection (48 h) cells were collected with 10 mM EDTA in 1xPBS and washed twice with cold PBS (samples used for phosphorylation analysis were collected and washed in the presence of phosphatase inhibitor cocktails). The cell pellets were stored at −80°C. Cells were lysed in lysis buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1 mM EDTA, 0.5% NP40, supplemented with protease and phosphatase inhibitor cocktails) for 30 min while rotating at 4°C and then centrifuge at high speed to collect the supernatant. The cell lysate was incubated with prewashed MagStrep beads (30 μl per reaction) for 3 h at 4°C. The beads were then washed four times with wash buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1 mM EDTA, 0.05% NP40, supplemented with protease and phosphatase inhibitor cocktails) and then proceed with on-bead digestion. The on-bead digestion was carried out as described before 39 . Briefly, after the final wash the beads were washed once with exchange buffer (50 mM Tris-HCl PH 7.5, 100 mM NaCl). Bead-bound proteins were alkylated with 3 mM iodoacetamide in the dark for 45 min and then quenched with 3 mM DTT for 15 min. Digestion was performed using 1 μl trypsin (1 μg/μl, Promega) overnight with shaking (1000 rpm) and the peptides were then purified. 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) (dilution used 1:1000). To purify clean 2xStrep-tagged N protein (mutant and control), we applied stringent washing and double elution strategy.
MS analysis using Orbitrap Fusion Lumos. The MS analysis was performed as described previously 82,83 with slight modifications. Briefly, approximately 0.5 µg of peptide mixture in 0.1% formic acid (FA) was injected into a nanotrap (PepMap 100, C18, 75 µm × 20 mm, 3 µm particle size) and desalted for 5 min with 0.1% FA in water at a flow rate of 5 µl/min. They were then eluted and analyzed using an Orbitrap Fusion mass spectrometer (MS) (Lumos, Thermo Fisher Scientific) coupled with an UltiMate™ 3000 UHPLC (Thermo Scientific). The peptides were separated by an EasySpray C18 column (50 cm × 75 µm ID, PepMap C18, 2 µm particles, 100 Å pore size, Thermo Scientific) with a 75-min gradient at constant 300 nL/min, at 40°C. The electrospray potential was set at 1.9 kV, and the ion transfer tube temperature was set at 270°C. A full MS scan with a mass range of 350-1500 m/z was acquired in the Orbitrap at a resolution of 60,000 (at 400 m/z) using the profile mode, a maximum ion accumulation time of 50 ms, and a target value of 2e5. The most intense ions that were above a 2e4 threshold and carried multiple positive charges (2-6) were selected for fragmentation (MS/MS) via higher energy collision dissociation (HCD) with normalized collision energy at 30% within the 2 s cycling time. The dynamic exclusion was 30 s. The MS2 was acquired with data type as centroid at a resolution of 30,000. Protein identification analysis from the raw mass spectrometry data was performed using the Maxquant software (v1.5.3.30) 38 as described 82 .
For phosphorylated peptides, we used Maxquant label-free quantification (LFQ) 38 . The analysis and quantification of phosphorylated peptides were performed according to published protocol 84 . The sample file naming format and the datasets associated with the MaxQuant analysis are represented in Supplementary Table S11.
Analysis of differential interaction. First, the identified protein group data were corrected for non-specific background binding by removing all proteins detected in mock control (cells transfected with the plasmid vector without N gene) affinity mass spectrometry (Supplementary Data 7). The normalized LFQ data were processed for statistical analysis on the LFQ-Analyst a web-based tool 85 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 cutoff of adjusted p-value < = 0.05, and Log foldchange > = 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 85 . BS 3 cross-linking. Bis(sulfosuccinimidyl) suberate (BS 3 , 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 86 . Briefly, 2xStrep-tagged N protein (mutant and control) was purified as mentioned above. The purified N protein (mutant and control) were cross-linked using 2 mM Bis(sulfosuccinimidyl) suberate (BS 3 ) for 30 min at room temperature. The control and cross-linked forms of N proteins (mutant and control) were separated on SDS-PAGE, subjected to silver staining and densitometry analysis of bands. GraphPad Prism (v9.1.1) was used for analysis and graph generation.
In vitro protein-RNA interaction (RIP) assay. The in vitro interaction assay was performed using purified 2xStrep-tagged N protein (mutant and control) and total isolated RNA from patient swabs as mentioned above in the RNA isolation section. The 2xStrep-tagged (in bead-bound condition) N protein (mutant and control) was incubated with total RNA in reaction buffer (50 mM Tris HCl pH7.5, 150 mM KCl, 0.1 mM EDTA, 1 mM DTT, 5% Glycerol, 0.02% NP40, 1 mM PMSF, 40 U RNase OUT™). After shaking incubation for 45 min, N proteins (mutant and control) were pulldown using MagStrep beads on a magnetic separator with four washes (using the same reaction buffer). After the final wash, the RNA was isolated by adding Trizol using the Zymo direct-zol method. The isolated RNAs were analyzed by RT-qPCR using specific viral N gene (N1 and N2), E gene, S gene, and ORF1ab region primers (Supplementary Table S10). GraphPad Prism (v9.1.1) was used for analysis and graph generation.
RNA-sequencing and differential gene expression analysis. Calu-3 cells were transfected with plasmids expressing the full-length N-control and N-mutant protein along with mock control. After 48-h cells were harvested in Trizol and total 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 Calu-3 RNA-sequencing were processed and trimmed using trimmomatic 63 and mapped to annotated ENSEMBL transcripts from the human genome (hg19) 87,88 using kallisto (v0.43.1) 89 . Differential expression analysis was performed after normalization using EdgeR integrated in the NetworkAnalyst 90 . GO biological process and pathway enrichment analyses on differentially expressed genes were performed using NetworkAnalyst 90 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Supplementary Information is available for this paper. Assembled virus genomes are available at GISAID (Supplementary Data 1). Raw viral reads, assembled virus genomes, and reads from RNA-Seq analysis of transfected Calu-3 cells have been uploaded to European Nucleotide Archive under the Study accession number PRJEB45515. The affinitypurification mass spectrometry (AP-MS) proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository 91 with the dataset identifier PXD027168. Source data are provided with this paper.