Increased flexibility of the SARS-CoV-2 RNA-binding site causes resistance to remdesivir

Mutations continue to accumulate within the SARS-CoV-2 genome, and the ongoing epidemic has shown no signs of ending. It is critical to predict problematic mutations that may arise in clinical environments and assess their properties in advance to quickly implement countermeasures against future variant infections. In this study, we identified mutations resistant to remdesivir, which is widely administered to SARS-CoV-2-infected patients, and discuss the cause of resistance. First, we simultaneously constructed eight recombinant viruses carrying the mutations detected in in vitro serial passages of SARS-CoV-2 in the presence of remdesivir. We confirmed that all the mutant viruses didn’t gain the virus production efficiency without remdesivir treatment. Time course analyses of cellular virus infections showed significantly higher infectious titers and infection rates in mutant viruses than wild type virus under treatment with remdesivir. Next, we developed a mathematical model in consideration of the changing dynamic of cells infected with mutant viruses with distinct propagation properties and defined that mutations detected in in vitro passages canceled the antiviral activities of remdesivir without raising virus production capacity. Finally, molecular dynamics simulations of the NSP12 protein of SARS-CoV-2 revealed that the molecular vibration around the RNA-binding site was increased by the introduction of mutations on NSP12. Taken together, we identified multiple mutations that affected the flexibility of the RNA binding site and decreased the antiviral activity of remdesivir. Our new insights will contribute to developing further antiviral measures against SARS-CoV-2 infection.


Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) was first discovered in 2019 and quickly spread around the world [1]. Novel SARS-CoV-2 variants have since continued to emerge and the number of virus-infected cases repeats increases and decreases [2]. The clinical spectrum of SARS-CoV-2 infection ranges from mild to critical. While most infections present mild or minor symptoms (e.g. fever, cough, sore throat, malaise, headache, muscle pain, nausea, vomiting, diarrhea, loss of taste and smell), severe acute respiratory disease requires admission to intensive care [3][4][5]. The illness can be observed even after successful vaccination [6]. Antiviral drugs that can be administered to patients after moderate or severe clinical symptoms have been observed have played important roles in clinical environments. Therefore, it is vital to understand the effectiveness of currently approved antivirals from multiple angles to develop future drugs. In particular, the potential to drive drug resistance should be evaluated because drug-resistant mutations have been observed in several viruses such as influenza A virus, human immunodeficiency virus and hepatitis B virus in the clinical environment [7][8][9][10].
Remdesivir (RDV) (GS-5734) is the US Food and Drug Administration (FDA)-approved drug for treatment of coronavirus disease 2019 (COVID-19) patients [11,12]. The compound is an intravenously administered adenosine analogue prodrug that binds to the viral RNAdependent RNA polymerase and inhibits viral replication. It has demonstrated antiviral activities against a broad range of RNA viruses including Ebolavirus, SARS-CoV, MERS-CoV and SARS-CoV-2 [13][14][15][16][17]. RDV has been widely used in the treatment of SARS-CoV-2 patients, however only two amino acid mutations (D484Y and E802D in non-structural protein [NSP] 12) were identified from SARS-CoV-2 patients that were administered RDV [18,19]. One mutation (E802D) was also found in in vitro serial passages of the virus under treatment of RDV [20]. Although studies regarding E802D revealed that the mutation decreased viral susceptibility to RDV [19,20], the mechanisms of how resistance arises have not yet been analyzed in detail. It is critical to elucidate the mechanisms of RDV resistance and to identify further RDV-resistant mutants that may arise in the future to circumvent resistance mutations before they become established in circulating strains.
To evaluate the effect of each gene mutation on viral propagation, genetically modified viruses should be engineered using the reverse genetics system. We recently established a quick reverse genetics system for SARS-CoV-2 using the circular polymerase extension reaction (CPER) method [21]. Nine viral genome fragments, which cover the full-length viral genome, and a linker fragment that encodes the promoter sequence were amplified by PCR and connected to obtain the circular viral DNAs by an additional PCR. By direct transfection of the circular DNAs, infectious SARS-CoV-2 was rescued. Introduction of reporters or mutations can be quickly completed by overlapping PCR or plasmid mutagenesis using the desired gene fragments of less than 5,000 base pairs (bp). While other reverse genetics systems for SARS-CoV-2 require specific techniques such as in vitro transcription or in vitro ligation, which are obstacles to mutagenesis [22,23], our method does not need these and has already been applied to the characterization of several viral mutations observed in the different SARS-CoV-2 variants [24,25], allowing us to simultaneously generate multiple mutants [26].
In this study, we attempted to identify multiple RDV-resistant mutations and examine the mechanisms of RDV resistance by a multidisciplinary approach that integrates state-of-the-art reverse genetics, mathematical modeling, and molecular dynamics analyses. We first predicted the presumed RDV-resistant mutations by in vitro passages of SARS-CoV-2 in the presence of RDV. Next, the recombinant viruses carrying the predicted mutations were generated by the CPER method and the efficiency of infectious virus production and antiviral effects of RDV on the mutants were examined by mathematical modeling. Finally, the conformational changes of NSP12 induced by mutations were analyzed by molecular dynamics simulations to understand the mechanisms of RDV resistance.

In vitro serial passages of SARS-CoV-2 in the presence of remdesivir
To identify the genes presumably involved in RDV resistance, we first passaged the SARS-CoV-2 strain SARS-CoV-2/Hu/DP/Kng/19-020 (GenBank: LC528232.2) in HEK293-C34 cells under treatment with RDV ( Fig 1A). In order to assess many viral mutations responsible for RDV resistance, we treated the cells with RDV at the highest concentration that would allow the virus to be passed on to the next generation. First of all, the cells were treated with 0.1 μM RDV, however cytopathic effect (CPE) was not observed in 14 days. The cells were then treated with 0.01 μM RDV 0.01 μM and CPE was observed at day 4. The culture supernatants were collected, stored as a Passage 1 (P1) sample ( Fig 1A). The concentration of RDV was gradually increased from 0.01 μM to 4.0 μM over 10 passages. Throughout the passages, the virusinfected cells were cultured until CPE was observed (3-8 days). Whereas no CPE was seen for 14 days after 0.1 μM RDV treatment in original viruses (P0), CPE was observed throughout the wells in the presence of 4.0 μM RDV at P10, indicating that the virus decreased susceptibility to RDV during the passages.
After 10 passages with RDV, the culture supernatants were collected (P10 with RDV) and subjected to MiSeq sequencing to determine the full-length viral sequence (DDBJ Accession number: LC742929). We prepared the viruses passaged 10 times without RDV as a control (P10 without RDV). Comparison of the P10 with RDV virus sequence with the original SARS-CoV-2 genome found 14 unique mutation sites and 6 mutation sites in P10 without RDV virus sequence ( Fig 1B and S1 Table). Five mutations found in P10 without RDV had a less than 50% mutant frequency, and mutations were not found in NSP7, 8, 12 nor 13, which are involved in the formation of the replication complex. Importantly, amino acid substitutions E796G and C799F in NSP12 were observed only in the P10 with RDV virus and there have been no reports of these two mutations to date, according to Nextstrain [27]. We then conducted the Sanger sequencing of the viruses after the limiting dilution cloning to investigate whether mutations are introduced in the same virus genomes or not. In total of 10 clones were obtained by the limiting dilution cloning after 10 times passage of SARS-CoV-2 with RDV (P10 with RDV). The deletion of nine nucleotides in NSP1 (82GHVM85V) was observed in 9 single viruses, and amino acid substitutions in NSP4 (V294L), NSP6 (L260F) and NSP12 (E796G and C799F) were observed in all 10 viruses, indicating that these mutations had been introduced in a single virus together (S2 Table). In addition, these mutations demonstrated more than 80% of frequency by Miseq sequence. According to Nextstrain, the same mutations had been detected in NSP1, NSP4 and NSP6, but only a few cases of each mutation had been reported. Besides, the deletion of 15 nucleotides in S, start-loss of ORF3A and S194P in N were observed in every single virus and detected by MiSeq sequence as well and these mutations have not been reported in Nextstrain.

Generation of RDV-resistant SARS-CoV-2
We then generated high-affinity NanoBiT (HiBiT)-carrying recombinant SARS-CoV-2 with each mutation to identify the RDV-resistant mutations based on strain JPN/TY/WK-521 (GISAID accession number: EPI_ISL_408667). NanoLuc enzymatic activity can be detected by interaction of HiBiT and large NanoBiT (LgBiT), which constitute a split reporter. The reporter SARS-CoV-2 can be generated by inserting only 11 amino acids into the viral genome, and HiBiT-carrying viruses exhibit similar growth kinetics to wildtype (WT) virus [21]. All recombinant SARS-CoV-2 with HiBiT and mutations were prepared using the CPER method that was previously established by our group. Amino acid substitutions were introduced by overlapping PCR and the full-length sequences of the mutant viruses were confirmed prior to assay by Sanger sequencing.

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 Because RDV acts as a nucleoside analog and targets the RNA-dependent RNA polymerase (RdRp) of coronaviruses, including SARS-CoV-2, in the current study we focused on the mutations in NSP12. We generated recombinant viruses with the E796G or C799F mutations that were observed in our P10 serial virus passage (Table 1). We also prepared recombinant viruses with mutations which demonstrated more than 80% of frequency in MiSeq sequence and was common in more than 9 out of 10 single virus clones after P10 (R10/E796G/C799F) or with mutations except for E796G (R10/C799F). In addition to the mutations observed in this study, we also characterized mutations that have been reported as, or anticipated to be, resistant to RDV, as listed in Table 1.
The amino acid mutation E802D in NSP12 was found during the serial passage of SARS-CoV-2 in vitro in the presence of RDV and another report showed that the same mutation was found in patients receiving RDV [19,20]. The D484Y mutation was also identified in a COVID-19 patient receiving RDV treatment [18]. Previously, amino acid substitutions F476L and V553L were identified as RDV-resistant mutations in the Betacoronavirus murine hepatitis virus (MHV) (15). The two affected amino acid residues (476F and 553L in MHV) are conserved across coronaviruses and correspond to 480F and 557V in the SARS-CoV-2 genome. We attempted to generate recombinant SARS-CoV-2 with either or both mutations but the virus with the single V557L mutation could not be rescued.
We then examined the sensitivity of recombinant viruses to RDV (S1 Fig). Viruses were cultured in the presence of RDV at 0-1.0 μM final concentration for 48 hours and luciferase activity was measured and normalized against control without RDV treatment (0 μM final concentration). The 50% effective concentration (EC 50 ) was calculated using the drc package (v3.0-1; R Project for Statistical Computing). All tested mutant viruses showed greater EC 50 than WT virus, although the difference between WT and F480L mutant was small, indicating that the mutations observed in NSP12 led to decreased viral sensitivity to RDV.

Time course analyses of infection with presumed RDV-resistant mutants
To investigate the growth efficiency of mutant viruses and the antiviral effects of RDV, we first performed time course analyses of infectious virus production with or without RDV for 96 hours and demonstrated the data up to 72 hours post infection (hpi) because most of the cells have been died and the viral titer have been decreased (Fig 2A and 2B). At all the indicated time points, the infectious titer of each mutant virus was similar or lower than that of the WT

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 virus in the absence of RDV treatment, indicating that the mutant viruses produced the infectious viruses with the same or lower efficiency as WT virus (Fig 2A). Conversely, significant differences were observed in the infectious titers of mutant viruses in the presence of RDV at 0.05 μM final concentration. In the left panel of Fig 2B, the infectious titers of mutant viruses (E796G, C799F, R10/E796G/C799F, and R10/C799F) gradually increased for 48 hours and were significantly higher than that of WT virus at 72 hpi. In the right panel of Fig 2B, the infectious titer of the E802D mutant virus increased rapidly and was significantly higher than that of WT virus at 48 and 72 hpi. The titers of the F484Y and F480L/V557L mutants were also significantly higher than that of WT virus. Meanwhile there were no differences between the titers of F480L mutant and WT viruses at the indicated time points. These results suggest that the sensitivity of all the mutant viruses, except for the F480L virus, to RDV was diminished, which was consistent with the results of the RDV susceptibility test demonstrating minimal change in the EC 50 of the F480L virus (S1 Fig).
Next, we investigated the the ratio of the virus-infected cells (Fig 2C and 2D). HEK293-C34 cells were infected with mutant viruses with and without RDV treatment. Virus-infected cells were harvested and fixed from 12-96 hpi and subjected to immunofluorescent assay using anti-SARS-CoV-2 NP antibody and DAPI. The virus infection rates were then calculated and demonstrated up to 72 hpi. All the mutant viruses demonstrated equivalent or significantly lower virus infection rates compared with WT virus in the absence of RDV. These data suggested that the number of cells infected with mutant viruses increased more slowly compared with WT virus in the absence of RDV treatment, consistent with the data on production of infectious virus particles. Meanwhile, the infection rates of the presumed RDV-resistant mutant viruses, except for D484Y, were higher or significantly higher (R10/E796G/C799F at 48 and 72 hpi, and R10/C799F at 72 hpi) than those of WT virus at 48 and 72 hpi in the presence of RDV, indicating that these mutant viruses can spread more efficiently than WT virus in the presence of RDV.

Antiviral effect of RDV on the presumed RDV-resistant mutants, analyzed by mathematical modeling
To quantify the kinetic parameters of SARS-CoV-2 and the antiviral effect of RDV on WT and RDV-resistant viruses, we developed a mathematical model for SARS-CoV-2 infection under RDV treatment. We examined the growth rate of HEK293-C34 cells up to 48 hours after seeding (S2A Fig), the degradation rate of SARS-CoV-2 at 37˚C (S2B Fig), the infectious virus production rate for 96 hours (Fig 2A and 2B), and the rate of infection in susceptible cells (Fig 2C  and 2D). These estimated parameters were fixed and used here.
To consider the variability of kinetic parameters and model predictions, we performed Bayesian estimation for the whole dataset using Markov chain Monte Carlo (MCMC) sampling, and simultaneously fit Eqs (1-4) with RDV (ε>0) and without RDV (ε = 0) to the concentrations of target cells, infected cells, and infectious virus (see Materials and Methods and S3 Fig). The estimated parameters are listed in Tables 2 and 3, and the simulation results of the model using these best-fit parameter estimates are shown with the data in S3 Fig. Comparing the virus production rate and the antiviral effect for the eight different resistant mutants (D484Y, F480L, F480L/V557L, E796G, C799F, R10/C799F, R10/E796G/C799F, E802D), we found that the virus production for all mutants was lower than that of the WT virus except for D484Y mutant (Fig 3B and Table 3). The D484Y mutation had a competitive advantage in virus production rate, and this property might be involved in its decreased susceptibility to RDV, although the difference from WT was only 1-1.25-fold.

Remdesivir resistance of SARS-CoV-2
Interestingly, all the other tested mutations strongly suppressed virus production. RDV showed more than 70% antiviral effect on three mutations, 81% (95% CI: 69-89) for D484Y; 80% (95% CI: 70-87) for F480L; and 78% (95% CI: 71-86) for F480L/V557L. However, the antiviral effect on two mutations was less than 40%, at 38% (95% CI: 29-49) for E802D and 38% (95% CI: 27-54) for R10/E796G/C799F ( Fig 3C and Table 3). All the examined mutations reduced the antiviral effect of RDV and the change was greatest in the mutations found in in vitro passages. The antiviral effect of RDV in the F480L mutation was not much different from that observed with WT virus, which is consistent with the results of the RDV susceptibility test. Because there was no correlation between the efficiency of virus production and the antiviral effect, the mechanisms of RDV resistance were predicted to be a structural defect in the direct interaction between the viral genome replication complex and RDV.

RDV-resistance mechanisms of the presumed RDV-resistant mutants, analyzed by computer simulation
Molecular modeling and molecular dynamics simulations were performed to clarify the structural and property changes caused by amino acid mutations in the NSP12 protein. The representative complete structure of the prepared protein-RNA complex is shown in Fig 4A. In this figure, NSP12, binding RNA, and RDV are shown in cartoon, ball and stick, and van der Waals notation, respectively. RDV is located at the end of the binding RNA and is inside the protein. Then, the root mean square deviations (RMSDs) were compared using molecular dynamics simulation trajectories of each complex of WT or mutant NSP12 protein and RNAincorporated RDV, as shown in Fig 4B. Most complexes reached thermodynamic stability and

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 plateau in RMSD after 200,000 steps. However, only the V557L mutant structure failed to reach the stabilized structure. When RNA structures for this mutant were superimposed and RMSDs were calculated for RNA and protein separately, the increase in RMSDs for RNA reached a plateau, while the RMSDs for protein continued to increase as before. This behavior means the bonds between protein and RNA tended to move apart. In other words, the RNAprotein complex tended to be unstable, which may correspond with the inability of this mutant virus to multiply, indicating the probable reason for failed rescue of the V557L mutant recombinant SARS-CoV-2 (Table 1). Therefore, we proceeded with the analysis of the other mutants only.
Using the trajectories obtained from molecular dynamics calculations, we calculated and compared the variation of RMSD for each substructure. The RMSDs of each mutant complex and the WT complex were searched for areas where they differed significantly, and the results are shown in Fig 4C. Locations where the RMSD variation of the mutant is greater than that of

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 the WT are indicated by pink spheres, and conversely, locations where the RMSD variation of the WT is greater than that of the mutant are indicated by blue spheres. In the tested mutations, except for F480L, the molecular vibration of the mutants tended to increase around the RNA-binding site as shown in S4 Fig, indicating that introduction of the mutations increased the flexibility of the RNA binding site. However, in the center of the RNA-binding site (near the RDV-binding site) of the F480L mutant, the molecular vibration of the mutant tended to be small, which is consistent with the small change in antiviral effect ( Fig 3C) and RDV sensitivity (S1 Fig). Taken together, NSP12 mutations found in previous studies and in our in vitro virus passages decreased the antiviral effect of RDV, though not to the same degree, and influenced increased flexibility of the RNA-binding site of the NSP12 structure.

Syrian hamster
Finally, we also examined the pathogenicity and viral growth kinetics of WT virus and virus with RDV-resistant mutations in in vivo hamster (Fig 5). The hamster group inoculated with WT virus showed significant weight loss compared to hamster inoculated with R10/E796G/ C799F virus at 4-7 dpi (Fig 5B). The oral swab samples were collected, and the viral RNA copies were compared for 5 days. In each time point, the viral RNA copies were significantly lower in the hamster group inoculated with R10/E796G/C799F virus. These findings indicate that R10/E796G/C799F virus has a less viral propagation efficiency in hamster than WT virus, which is consistent with the lower virus production rate of R10/E796G/C799F virus in HEK293-C34 cells (Figs 2A, 2C and 3B).

Discussion
As demonstrated by the recent emergence of the Omicron strain, mutations have continuously accumulated on the SARS-CoV-2 genome [28]. To implement the most effective measures against COVID-19, it is crucial to predict clinically important mutations in advance, including drug resistance mutations. In this study, we generated eight recombinant viruses with presumed RDV-resistant mutations, using a simple reverse genetics system. Quantifying the kinetic parameters for RDV-resistant viruses showed no dramatic increases in the efficiency of infectious virus production in any mutant virus. Besides, R10/E796G/C799F virus was less efficient in virus propagation in hamsters than WT virus. But importantly, all mutants rendered RDV ineffective to some extent. Molecular modeling and molecular dynamics simulations of the mutant NSP12 proteins revealed that the tested mutations, excluding F480L, contributed to increased molecular vibrations around the RNA-binding site. Our multidisciplinary approach of molecular virology, mathematical, and molecular modeling discovered mutations involved in RDV resistance and the mechanisms of this drug resistance.
When evaluating virus growth kinetics, titers of different viruses are generally compared at the same time points. However, the total numbers of cells susceptible to viral infection are not constant between wells. The faster a virus propagates and spreads, the faster the number of cells available for virus infection decreases. Thus, comparison of titers at specific time points may not be the best method of evaluating the proliferative capability of viruses with different properties or the effects of drugs against them. Here we combined mathematical models and statistical methods, analyzed intercellular infection dynamics of SARS-CoV-2 in the momently changing cells, and identified the differences in infectious virus production and antiviral effects between WT and mutant viruses (Fig 3).

Remdesivir resistance of SARS-CoV-2
Furthermore, molecular simulation of the NSP12 structure clarified the mechanism causing failure of V557L to proliferate and provided new insights into the mechanisms of RDV resistance of the tested mutants (Fig 4), in contrast to other previous studies that only identified the  S5C Fig). Besides, RMSD for the average structure of the whole model after 300,000 steps, which is considered to be equilibrated, was 1.64 Å. These findings indicates that thermodynamic oscillations of NSP12 is not changed largely, and the small model which was used in this study is accurate enough to evaluate the effects of the mutations on molecular vibration. In mutants resistant to RDV, regardless of the location of the mutation, large differences were observed in the RNA binding site when thermodynamic oscillations were compared with WT. Mutations not involved in resistance showed the opposite trend of variation in terms of protein flexibility. We therefore speculate that flexibility in the binding site may be a factor in resistance to RDV. More detailed mechanisms of resistance might be revealed by evaluating polymerase activities of mutants.
Three mutations (C799F, E796G, and E802D) markedly diminished the antiviral effect of RDV (Fig 3C). These mutations did not contribute to efficient virus production but increased the flexibility of the RNA-binding site, which probably enabled these viruses to evade the functional inhibition by RDV binding. The C799F virus was less efficient in virus production than R10/C799F virus ( Fig 3B). Besides, the mutations introduced in R10/C799F virus was not detected in the viruses passaged 10 times in absence of RDV by MiSeq sequeincig and observed in the same individual virus genomes (S1 and S2 Tables). Therefore, some mutations, which were introduced in R10/C799F or R10/ E796G/C799F and observed in genome regions other than NSP12, might restore the efficiency of virus production in the mutant viruses. Unlike the mutations found in in vitro analyses, D484Y, F480L, and F480L/V557L affected neither virus production nor RDV sensitivity. We expect that even the homologous amino acid mutations (F480L and V537L) would have different effects on the NSP12 of MHV and SARS-CoV-2, since there are also many different amino acid residues within NSP12 that interact with them to change the structure of NSP12. Importantly, these results consistent between in vitro passages and computer simulations highlight the accuracy and usefulness of current interdisciplinary research.
In the previous study, 2 amino acid substitutions were observed in NSP12 by in vitro passages of MHV with remdesivir and the mutations decreased viral fitness of MHV in vitro. Mutant SARS-CoVs with the corresponding mutations also attenuated the pathogenesis in mice [15]. Therefore, we anticipated that the NSP12 mutations observed in in vitro passage of SARS-CoV-2 with the presence of remdesivir in this study had neither gained pathogenicity nor fitness and then constructed them by CPER method. Every mutation observed in in vitro passages of SARS-CoV-2 failed to increase the efficiency of infectious virus production. Morevover, R10/E796G/C799F virus was less effective than WT virus in propagation in hamsters. Co-infection competition assay in previous studies revealed that E802D on SARS-CoV-2 or F476L/V553L mutations on MHV decreased fitness [15,20], consistent with our results. Although further characterization is necessary, these results indicate that the observed mutations in this study are unlikely to be persistent in the virus population without RDV or dramatically accelerate the SARS-CoV-2 epidemic. However, because the E802D mutation found initially during in vitro passage was actually detected in clinical patients [19,20], it is quite possible that the analyzed mutations may be reported following the sustained administration of RDV, and thus continuous viral sequence analyses from SARS-CoV-2 patients treated with RDV is vital. Currently, we cannot definitively state that the repeated administration of RDV will be problematic, but additional compounds with higher affinity for the RdRp complex than RDV are likely to become desirable.
The SARS-CoV-2 pandemic is not coming to an end, rather, new virus strains have sequentially emerged. Deletions and mutations that can facilitate higher transmissibility or antibody

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 evasion have been highlighted in the Omicron variant [29][30][31][32][33]. It has been revealed that two doses of mRNA vaccine are insufficient to provide immunity against the Omicron variant [34][35][36], and new vaccines that will provide more robust immunity are urgently needed. In addition, new antivirals have been developed and are expected to be approved for clinical use [37][38][39]. While development of new vaccines and antivirals will continue, it remains important to evaluate their safety, and to pay attention to resistant mutations. Our state-of-the-art study, in which the drug sensitivities of multiple mutant viruses were simultaneously determined, will accelerate the development of new measures, evaluation of drug resistance and deepen our understanding of the driving forces for mutation of SARS-CoV-2.
SARS-CoV-2 strain SARS-CoV-2/Hu/DP/Kng/19-020, was kindly provided by Dr. Sakuragi at the Kanagawa Prefectural Institute of Public Health and strain JPN/TY/WK-521 by Dr. Masayuki Shimojima at the National Institute of Infectious Diseases. All experiments involving SARS-CoV-2 were performed in biosafety level 3 laboratories, following the standard biosafety protocols approved by the Research Institute for Microbial Diseases at Osaka University.

Chemical inhibitors and antibodies
RDV (GS-5734) was purchased from Cayman Chemical, dissolved in dimethyl sulfoxide (DMSO) and stored at 50 mM at −30˚C. To detect SARS-CoV-2-infected cells, mouse monoclonal antibody against SARS-CoV-2 NP (Clone# S2N4-1242) was kindly provided by Bio Matrix research. Alexa Fluor 488-conjugated anti-mouse antibodies were purchased from Life Technologies.

Serial passages of SARS-CoV-2
SARS-CoV-2 strain SARS-CoV-2/Hu/DP/Kng/19-020 was serially passaged in HEK293-C34 cells 10 times in the presence of RDV (Fig 1). HEK293-C34 cells were prepared with DMEM containing 10% FBS, blasticidin, and 1 μg/ml doxycycline hydrochloride in six-well plates. One day later, the cells were infected with SARS-CoV-2 at a multiplicity of infection (MOI) of 0.01 for 1 hour to allow virus attachment. Culture supernatants were then replaced with DMEM containing 2% FBS, blasticidin, 1 μg/ml doxycycline hydrochloride, and RDV. The virus-infected cells were incubated and the supernatants were collected when CPE was observed throughout the wells. The collected supernatants were centrifuged at 1,500 ×g for 5 min to remove cells and debris, and 10 μl of the supernatants were passaged. The titer of the virus in each passage was demonstrated in S4 Table. The final concentration of RDV was gradually increased from 0.01 μM in P1 to 4.0 μM in P10. As a control, SARS-CoV-2 strain SARS-CoV-2/Hu/DP/Kng/19-020 was also passaged in HEK293-C34 cells 10 times in the absence of

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 RDV. The virus was infected at a MOI of 0.01 and the culture supernatants were collected at day 2. The collected supernatants were centrifuged at 1,500 ×g for 5 min to remove cells and debris, and 10 μl of the supernatants were passaged.

Validation of the virus genome sequence using Sanger sequence
Total RNA was extracted from the supernatants of SARS-CoV-2-infected cells by using a Pure-Link RNA mini kit (Invitrogen) and subjected to cDNA synthesis using a PrimeScript RT reagent kit (Perfect Real Time) (TaKaRa Bio) and random hexamer primers. A total of nine DNA fragments, covering the full-length SARS-CoV-2 genome were amplified by PCR using a PrimeSTAR GXL DNA polymerase (TaKaRa Bio), the synthesized cDNA and specific primer sets from CoV-2-G1-Fw to CoV-2-G10-Rv designed previously [21]. The amplified PCR fragments were purified using a gel/PCR DNA isolation system (Viogene) and sequenced in both directions using the ABI PRISM 3130 genetic analyzer (Applied Biosystems) with specific primers for SARS-CoV-2.

Limiting dilution cloning
After the serial passage of SARS-CoV-2 10 times with RDV (P10 with RDV), the virus was prepared at 10 3 TCID 50 /ml. Then virus was conducted three-fold serial dilution with DMEM containing 2% FBS and inoculated into VeroE6/TMPRSS2 cells prepared in 96-well plates. After 72 hours, the supernatants were collected from a well where a single focus was observed. The collected supernatants were centrifuged at 3,000 ×g for 5 min to remove cells and debris and the virus genome was validated by the Sanger sequence.

Rescue of the presumed RDV-resistant viruses
All the HiBiT-carrying SARS-CoV-2 with RDV-resistant mutations were rescued by the CPER method, which was established in our previous study [21]. Briefly, nine cDNA fragments, covering the entire genome of SARS-CoV-2 were prepared by PCR using PrimeSTAR GXL DNA polymerase and SARS-CoV-2 viral gene fragment-encoding plasmids. In addition, an untranslated region (UTR) linker fragment encoding the 3 0 43 nucleotides (nt) of SARS-CoV-2, hepatitis delta virus ribozyme (HDVr), bovine growth hormone (BGH) poly(A) signal, cytomegalovirus (CMV) promoter, and the 5 0 25 nt of SARS-CoV-2, was amplified by PCR. A HiBiT luciferase gene (VSGWRLFKKIS) and a linker sequence (GSSG) were introduced into the N terminus of the ORF6 gene of SARS-CoV-2 by site-directed mutagenesis of the viral genome fragment-cloning plasmid and the plasmid was used as a template to amplify a cDNA fragment. Presumed RDV-resistant mutations were introduced into the SARS-CoV-2 cDNA fragments by overlap PCR using specific overlapping primer sets (S3 Table). The nine

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 SARS-CoV-2 cDNA fragments and the UTR linker fragment (0.1 pmol each) were mixed together and subjected to CPER. The CPER products were then directly transfected into HEK293-C34 cells using Trans IT LT-1 (Mirus). At 6 hours post-transfection, the culture media were changed to DMEM containing 2% FBS, blasticidin, and doxycycline hydrochloride (1 mg/ml). When CPE was observed throughout the wells (usually around 7 days post-transfection), the culture supernatants were collected and centrifuged at 1,500 ×g for 5 min to remove cells and debris. The culture supernatants were then passaged once using VeroE6/ TMPRSS2. The virus sequences were confirmed in the passaged virus solutions by Sanger sequencing and thereafter the virus stocks were stored at −80˚C until use. The experimental procedure for the construction of mutant SARS-CoV-2 was approved by the committee on genetically modified organisms in Osaka university (project number 4680) and Japanese Minister of Education, Culture, Sports, Science and Technology in compliance with the Cartagena Protocol on Biosafety (2受文科振第764号). All experiments using mutant SARS-CoV-2 were performed in fully licensed BSL-3 facilities in Osaka university.

Virus titration
Infectious titers in culture supernatants were determined by 50% tissue culture infective doses (TCID 50 ). The TCID 50 was calculated with Reed-Muench Method [46]. VeroE6/TMPRSS2 cells were prepared in 96-well plates and infected with SARS-CoV-2 after ten-fold serial dilution with DMEM containing 2% FBS. Virus titers were determined at 72 hpi.

RDV susceptibility analysis using the HiBiT system
HEK293-C34 cells were seeded in 48-well plates in DMEM with 10% FBS, blasticidin and 1 μg/ml doxycycline hydrochloride. One day later, HiBiT-carrying SARS-CoV-2 was allowed to attach for 1 hour. Culture media were then replaced with new media containing 2% FBS, blasticidin, 1 μg/ml doxycycline hydrochloride, and RDV (0-1.0-μM final concentration). At 48 hpi, luciferase activity was measured using a Nano-Glo HiBiT lytic assay system (Promega), following the manufacturer's protocols. Briefly, Nano-Glo substrate including LgBiT protein was added to the virus-infected cell lysates after all culture supernatants were removed. Luciferase activities were measured using a luminometer and normalized to luminescence without RDV treatment (0-μM final concentration). The EC 50 was calculated using the drc package (v3.0-1; R Project for Statistical Computing).

Time course analyses of infectious virus production
HEK293-C34 cells were prepared in 96-well plates in media containing 1 μg/ml doxycycline hydrochloride. Cells were infected with HiBiT-carrying viruses at MOI = 0.01 for 1 hour. Culture media were changed to fresh media containing 2% FBS, blasticidin, and 1 μg/ml doxycycline hydrochloride, with or without 0.05 μM RDV (final concentration). At 12,24,48,72, and 96 hpi, culture supernatants of the virus-infected cells were collected and infectious titers in the supernatants (TCID 50 /ml) were determined by virus titration.

Time course analyses of infection rates in cells
After removal of the culture supernatants at the indicated time points in the time course analyses of infectious virus production, the cells were fixed with 4% paraformaldehyde (Nacalai Tesque). The fixed cells were permeabilized with 0.2% Triton X-100 (Nacalai Tesque) in PBS for 20 min, blocked with 1% bovine serum albumin fraction V (Sigma) in PBS, and then reacted with anti-SARS-CoV-2 NP antibody in PBS for 1 hour at room temperature. After washing

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 with PBS three times, the cells were incubated with a 1:1,000 dilution of goat anti-mouse IgG Alexa Fluor 488-conjugated secondary antibody (Thermo Fisher Scientific) in PBS for 1 hour at room temperature. The cells were then incubated with DAPI (Thermo Fisher Scientific) (1:2,000 dilution) for 10 min. Immunopositive signals were confirmed under a FluoView FV1000 confocal laser scanning microscope (Olympus), with appropriate barrier and excitation filters. Quantitative imaging data were obtained using a CellVoyager CQ1 benchtop highcontent analysis system (Yokogawa Electric Corporation) and analyzed with CellPathfinder high content analysis software (Yokogawa Electric Corporation). The number of SARS-CoV-2-infected cells stained by anti-SARS-CoV-2 NP antibody and the number of cell nuclei stained by DAPI were counted. The infection rates were then calculated by dividing the number of SARS-CoV-2-positive cells by the total number of cell nuclei.

Growth of HEK293-C34 cells
HEK293-C34 cells were seeded in 48-well plates with DMEM containing 10% FBS and blasticidin. At 24 hours post-seeding, media were replaced with DMEM containing 10% FBS, blasticidin, and 1 μg/ml doxycycline. All the cells were collected, and the cell numbers were counted every 12 hours for 48 hours post-medium change.

Degradation rate of HiBiT-carrying SARS-CoV-2
HiBiT-carrying SARS-CoV-2 were incubated at 37˚C with 5% CO 2 for 48 hours. Every 12 hours, the virus solutions were collected and subjected to virus titration to quantify infectious virus.

Quantification of cell growth and virus decay kinetics
To estimate the growth kinetics of target cells, we used the following mathematical model: where the variable T(t) represents the number of uninfected target cells (cells/ml) at time t, and the parameters g and K indicate the growth rate and the carrying capacity of the target cells (cells/ml), respectively. Using the non-linear least square method, we fitted the model to the time-course growth data of cells (see Growth of HEK293-C34 cells and in S2A Fig) and estimated g and K.
Furthermore, we estimated the clearance rate of infectious viruses, c, by a simple exponential decay model: where V(t) represents the amount of infectious virus (TCID 50 Table 2.

Mathematical model for SARS-CoV-2 infection
We employed the following mathematical model for SARS-CoV-2 infection in cell culture considering the antiviral efficacy of RDV: where T(t), E(t), and I(t) are the numbers of uninfected target cells, eclipse phase cells, and virus-producing cells (cells/ml) at time t, respectively, and V(t) is the amount of infectious virus (TCID 50 /ml) at time t. The uninfected target and eclipse phase cells divide in logistic manner at rate g and carrying capacity K. The target cells are infected by viruses at rate β, and the virus-infected cells stay in the eclipse phase during the period 1/k. After this, they become virus-producing cells. The progeny viruses are produced by the virus-producing cells at rate p.
The parameters δ and c indicate the death rate of infected cells and the clearance rate of viruses, respectively. The inhibition rate of virus production by RDV is assumed to be ε. The fold-change of virus production rates of RDV-resistant viruses compared with WT virus are η (i.e., η = 1 for WT virus).

Data fitting and parameter estimation
The parameters g, K, and c were independently estimated and fixed. A statistical model adopted in Bayesian inference assumed that measurement error followed a normal distribution with mean zero and constant variance (error variance). A gamma distribution was used as a prior distribution, and it inferred a distribution of error variance. As an output of MCMC computations, the posterior predictive parameter distribution represented parameter variability, and it inferred distributions of model parameters and initial values of variables. The estimated parameters and initial values are listed in Tables 2 and 3.

Theoretical predictions and analyses of the effects of amino acid mutations
The WT NSP12 structure used as a reference was created based on the crystal structure (PDB ID: 6XEZ) [47]. From this crystal structure, only the NSP12 protein and 28 bases of the binding RNA (binding site) were extracted. The target mutant structures were constructed by substituting amino acids in the above WT structure. Based on the crystal structure (PDB ID: 7BV2) [48], the terminal nucleotides of each predicted structure were replaced with RDV. To neutralize the charge of these complex structures, counter ions were placed, and sufficient water molecules were placed around them. Each structure was stabilized by the energy minimization method and used as the initial structure for molecular dynamics simulations. The composite structure was thermally stabilized by raising the temperature from 0 K to 310 K (in vivo temperature, approximately 36.85˚C) over 500,000 steps with Δt = 0.2 fs. The structural changes during this temperature increase process were structurally sampled for each complex at every 1000 steps. These sampling structures were superimposed, and the RMSD for each protein was calculated. In addition, to observe the extent to which the structural properties differed between WT and mutants, we superimposed them in various substructures and calculated and compared the RMSD differences. These energy minimizations and molecular dynamics simulations were performed with the AMBER18 program package [49]. The

Animal experiments
Syrian hamsters (male, 4 weeks old) were purchased from Japan SLC. Baseline body weights were measured before infection. For the virus infection experiments, hamsters were anesthetized by intramuscular injection of a mixture of 0.15 mg/kg medetomidine hydrochloride (Domitor, Nippon Zenyaku Kogyo), 2.5 mg/kg butorphanol (Vetorphale, Meiji Seika Pharma) and 4.0 mg/kg alphaxaone (Alfaxan, Jurox). WT and R10/E796G/C799F viruses (10 4 TCID 50 in 100 μl) or saline (100 μl) was intranasally inoculated under anesthesia. Oral swabs were collected at the indicated timepoints and body weight was recorded daily by 7 dpi. The viral RNA load in the oral swabs was determined by RT-qPCR. Total RNA was extracted from the oral swabs using a PureLink RNA mini kit (Invitrogen). The RNA was used as the template for RT-qPCR performed in accordance with the manufacturer's protocol using the One Step Prime-Script III RT-qPCR Mix (Takara) and the following primers and probe: Forward

Statistical analysis
Results are indicated as the means ± standard deviations or standard errors. Statistical significances were determined by two-tailed Student's t test, the one-way ANOVA with Dunnett's test or the Kruskal-Wallis test with the two stage linear step-up procedure of Benjamini, Kreiger, and Yekutieli, which was performed using GraphPad Prism (Software ver. 9.2.0). Significantly different values are indicated by asterisks (*p<0.05 or ***p<0.001). The RNA binding space is exposed by using

PLOS PATHOGENS
Remdesivir resistance of SARS-CoV-2 only the protein surface notation and cartoon notation.