Use of whole genome sequencing to identify low‐frequency mutations in SARS‐CoV‐2 patients treated with remdesivir

Abstract Background Remdesivir (RDV) has been shown to reduce hospitalization and mortality in COVID‐19 patients. Resistance mutations caused by RDV are rare and have been predominantly reported in patients who are on prolonged therapy and immunocompromised. We investigate the effects of RDV treatment on intra‐host SARS‐CoV‐2 diversity and low‐frequency mutations in moderately ill hospitalized COVID‐19 patients and compare them to patients without RDV treatment. Methods From March 2020 to April 2022, sequential collections of nasopharyngeal and mid‐turbinate swabs were obtained from 14 patients with and 30 patients without RDV treatment. Demographic and clinical data on all patients were reviewed. A total of 109 samples were sequenced and mutation analyses were performed. Results Previously reported drug resistant mutations in nsp12 were not identified during short courses of RDV therapy. In genes encoding and surrounding the replication complex (nsp6‐nsp14), low‐frequency minority variants were detected in 7/14 (50%) and 18/30 (60%) patients with and without RDV treatment, respectively. We did not detect significant differences in within‐host diversity and positive selection between the RDV‐treated and untreated groups. Conclusions Minimal intra‐host variability and stochastic low‐frequency variants detected in moderately ill patients suggests little selective pressure in patients receiving short courses of RDV. The barrier to RDV resistance is high in patients with moderate disease. Patients undergoing short regimens of RDV therapy should continue to be monitored.

The SARS-CoV-2 pandemic has resulted in millions of deaths, and while vaccines have been shown to have high effectiveness in preventing severe disease, 1 there have been breakthrough infections, and some populations remain at risk for severe disease.SARS-CoV-2 specific small molecule inhibitors offer a promising means of reducing hospitalizations and mortality.Remdesivir (RDV) is a nucleotide analog that is administered as a prodrug, and subsequently converted by the host cell into the active triphosphate form to inhibit viral replication 2 and has maintained activity against all variants that have emerged up to Omicron. 3 Therapeutic efficacy demonstrated by increased survival has been noted in hospitalized and out-patient trials. 4,5Resistance conferring mutations to RDV are very rare based on a survey of SARS-CoV-2 genomes available in GISAID 6 and has been primarily reported in immunocompromised patients who are on prolonged therapy. 7,8Additionally, a recent report highlighted an immunocompromised patient who had rebound infections and multiple courses of RDV, but no resistance mutations were identified. 9Overall, there is a paucity of data on patients with serial sample collections over the course of their treatment making it difficult to determine when, or if, resistance mutations arise, and at what frequency.
Inadequate therapy can be a driver for the selection of drug resistance.While mutations conferring resistance to RDV have been identified in patients who received prolonged therapies for SARS-CoV-2, 7,10 it remains unclear if hospitalized individuals who received shorter courses of therapy are likely to develop mutations that confer resistance to RDV.High-throughput sequencing is a method for detecting low-frequency variants that emerge following the administration of antiviral therapy. 11,12In this study, we collected serial samples from patients treated with RDV, including those who had shortened courses of therapy.Novel polymorphisms identified as lowfrequency variants were analyzed for their potential resistance profiles and tested for selection.

| Whole genome sequencing (WGS)
The ARTIC SARS-CoV-2 protocol (https://artic.network/ncov-2019)was used for amplicon generation.Sequencing was performed using the methods outlined. 13,14Briefly, nucleic acid extraction was performed 15 and cDNA was generated by the ARTIC protocol.After cDNA synthesis, two multiplex PCR tiling reactions were prepared using the ARTIC V3, V4, and V4.1 primers (depending on sampling date).cDNA was combined with Q5 High-Fidelity 2X Master Mix (NEB, USA).To mix #1, 10 μM of ARTIC primer pool #1 were added and 10 μM of ARTIC primer pool #2 were added to mix #2.Cycling was then performed as follows: 98 C for 30 s followed by 35 cycles of 98 C for 15 s and 65 C for 5 min.Following this reaction, products were combined, and a clean-up with AMPure XP beads was performed (Beckman Coulter, USA).The quantity of amplicons was measured with the Qubit 4.0 fluorometer using the 1X dsDNA HS Assay Kit (Thermo Fisher Scientific, USA).The sequencing libraries were prepared using the Nextera DNA Flex Prep kit (Illumina, USA) as per manufacturer's instructions.Paired-end (2 Â 150 bp) sequencing was performed on a MiniSeq with a 300-cycle reagent kit (Illumina, USA).

| Phylogenetic analysis
Transition (Ti) and transversion (Tv) substitutions, deletions, and insertions were tabulated from Breseq results.Hypothesis Testing using Phylogenies (HyPhy) v2.5.42 was used for testing signals of positive selection 24 in the ORF1ab gene on RDV-treated sequences.Prior to selection analysis, Virulign 25 v1.0.1 was used to generate codon alignments of all sequences with the SARS-CoV-2 Wuhan/Hu-1/2019 (MN908947.3)ORF1ab as the reference.A HyPhy inferred tree was generated and phylogenies were labeled using the Phylotree Widget.
Excess stop codons were removed using HyPhy before performing selection analysis.The adaptive branch-site random effects likelihood (aBSREL) 26 method from HyPhy was used to test for evidence of positive selection among the foreground RDV branches relative to the untreated sequences.

| Cohort description
In total, 44 patients who were hospitalized from March 2020 to April 2022 and had at least two SARS-CoV-2 positive nasal samples were included in the study.Patient demographic and clinical characteristics are shown in Table 1.Fourteen (31.8%) patients received RDV for a median duration of 4 days (IQR, 2.5-9).Patients were sampled a median of 2 (IQR, 2-3) times following commencement of therapy and collection times are shown in Table 2. Samples were subjected to WGS if SARS-CoV-2 RT-PCR cycle threshold values for the E-gene target were less than 30 (Tables S1 and S2

| Drug resistance mutations were not identified by WGS in RDV-treated patient cohort
Multiple mutations have been reported that confer varying levels of resistance to RDV in the viral replication complex. 6We examined 35 serial collections obtained from 14 RDV-treated patients.The consensus genomes were examined and neither E802D, a mutation in nsp12 associated with treatment failure, 7 nor other known resistance mutations were not detected as a consensus mutation in any patient samples regardless of the duration of RDV therapy.In parallel, we analyzed 74 serial samples collected from 30 patients who did not receive RDV therapy and did not identify any drug-resistance mutations.

| Low-frequency variants identified during treatment
The highly sensitive sequencing platform allows for detection and quantification of low-frequency variants that would not be identified  in the consensus sequence.Work by other groups suggests that only mutations present at >15% were likely to represent viral adaptation in the face of selective drug pressure, thus we used a similar cut off for classifying low-frequency mutations under this threshold. 29w-frequency minority variants within the replication complex and surrounding genes (nsp6-nsp14) were identified in 7/14 (50%) of RDV-treated patients (Figure 1A).These low-frequency variants were unique instances and not repeatedly detected across patients.Minority variants were identified in 18/30 (60%) of patients without RDV treatment (Figure 1B).Low-frequency variants identified in early time points did not persist in subsequent collections and did not become fixed in either group.Additionally, no low-frequency variants were observed at E802 or V166, both positions associated with RDV resistance in patients. 7,10Low-frequency variants were also observed in other regions of the genome in both treatment groups (Figure S1a and b).

| Within-host variability detected in RDVtreated and untreated patient cohorts
We investigated intra-host viral diversity for patients treated with and without RDV therapy.Changes in the total number of amino acid mutations identified by Breseq, including synonymous and nonsynonymous single nucleotide polymorphisms (SNPs), insertions, deletions, and complex substitutions across samples from different time points can reveal intra-host viral diversity.Intra-host variability was detected in nsp6-nsp14 across serial collections in 10/14 (71.4%) patients treated with RDV (Figure 1A).However, there were no changes in the consensus sequences of serial collections in 4/14 (28.6%) patients, indicating no intra-host diversity.Within-host diversity in nsp6-nsp14 was detected in 20/30 (66.7%) patients without RDV therapy (Figure 1B).There were 13/14 (92.8%) and 27/30 (90.0%) patients in the RDV-treated (Figure S1a) and untreated (Figure S1b) groups, respectively, where intra-host variation across serial collections were detected throughout the whole genome.Total within-host diversity was also compared among the RDV-treated and untreated groups by investigating the absolute difference in the number of total mutations per sample relative to the first analyzed serial collection (Figure 2).
Within-host diversity between the RDV-treated and untreated groups in the nsp6-nsp14 region or the whole genome were not significantly different based on the Mann Whitney U test.
3.5 | Within-host variability and selection processes are minimal in RDV-treated and untreated patients Transition (Ti) and transversion (Tv) events were examined across the whole genome for patients treated with and without RDV therapy.
Transitions are substitutions interchanging bases of the same ring structure while transversions are interchanges of purine bases for pyrimidine bases.The total number of all transition and transversion events from each sample from all patients among the RDV-treated and untreated group were compared for the nsp6-nsp14 region (Figure 3A) and the whole genome (Figure 3B).C ! T Ti were the most common type of substitutions in both groups.Prior to correcting for phylogenies, significant differences in the number C ! T mutations in the nsp6-nsp14 region (Figure 3A) and A ! G, C ! T, G ! A, G ! T, and T !C (Figure 3B) mutations in the whole genome between the RDV-treated and untreated group were observed.However, these differences can be attributed to the varying makeup of lineages in each treatment group, as the untreated group contains mostly early lineages.In the RDV-treated group, 1165 Ti and 805 Tv were identified within the whole genome, giving a Ti/Tv ratio of 1.45.
There were 1245 Ti and 664 Tv in the untreated group within the whole genome with a Ti/Tv ratio of 1.88.When comparing the Ti/Tv ratios among both groups, no significant difference was observed.
Furthermore, we identified 144 small and mid-size deletions across the whole genome in the 14 RDV-treated patients of which 18 were low frequency (Figure S1a).142 deletions were detected in the 30 patients without RDV treatment with 43 low-frequency deletions Absolute difference in number of total consensus and low-frequency mutations (including SNPs, amino acid insertions, amino acid deletions, and complex mutations) relative to the first examined sample of each patient within nsp6-nsp14 and the whole genome.The Mann Whitney U test was used to calculate p-values comparing distributions between RDV-treated and untreated groups.
(Figure S1b).Deletions are denoted on heatmap figures by "Δ bp coding."Most deletions detected were found in the consensus sequences of samples and are characteristic of VOCs particularly the Delta and Omicron lineages.
To determine whether selective pressures were exerted in the RDV-treated group, we tested for episodic positive selection while correcting for the phylogenetic structure of our samples.Codon alignments for ORF1ab were generated for all sequences and the aBSREL model 26 from HyPhy 24 was applied to an RDV-treatment labeled phylogeny generated by the Phylotree Widget.No evidence of episodic diversifying positive selection was identified in the ORF1ab gene in the foreground RDV branches relative to the untreated branches.
Together, this suggests that evolutionary processes such as positive selection, replicase complex adaptation, and resistance emergence are minimal with RDV treatment.

| DISCUSSION
In this study, WGS was used to identify low-frequency mutations in patients with serial collections during RDV therapy.7][8][9] Overall, RDV did not result in the generation of low-frequency drug-resistance mutations in the are advantageous to the virus may exist at low-frequency and emerge under selective pressure. 10Our data suggest that if this is the case, strong selection is not present in the early phase of treatment.We speculate that minimal intra-host diversity and the emergence of lowfrequency variants arise stochastically due to genetic drift and not from RDV selective pressures.Other studies in persistently infected immunocompromised patients have also failed to identify resistance mutations to RDV, 30 although interestingly, mutations that reduce monoclonal antibody neutralization efficacy arise more easily.
Our study is novel as it focused on patients who received less than 10 days of therapy and indicates that this does not select for resistant viruses in hospitalized patients.While it cannot be discounted that mutations arose at later time points, this seems unlikely as all RDV-treated patients recovered and none had prolonged infection.This is indirectly supported by the fact that RDV has been used as 3-day treatment for non-hospitalized patients and has not resulted in documented resistance emerging. 5Collectively, these findings suggest that the barrier to resistance is high and that hostimmunosuppression is needed for the development of resistance.We highlight these effects in moderately ill patients with shorter-term infection without severe immunocompromise.This underscores the relative safety of RDV against the emergence of drug resistant mutations in such patients.Our study had several limitations, including the small cohort size, and that only patient samples with sufficiently high viral load (e.g., Ct < 30) were suitable for sequencing.The phylogenetic makeup of our treatment groups also varied.We did attempt to control for the phylogenetic structure while testing for positive selection; however, this will still reduce our statistical power.
Here, we highlight the effects of short courses of RDV on lowfrequency mutations and SARS-CoV-2 intra-host variability in patients with moderate disease.Resistance mutations were not detected within the replication complex and intra-host diversity remains minimal in RDV-treated patients.This suggests that the barrier to resistance is high, likely requiring prolonged infection and immunosuppression of the host.Future improvements to wholegenome sequencing assays are needed to interrogate samples with low viral loads and improve surveillance for antiviral resistance mutations.
Serial samples obtained from hospitalized adult patients with laboratory-confirmed COVID-19 were recruited from five hospital sites in the greater Toronto area.Patients were prospectively recruited by the Toronto Invasive Bacterial Disease Network (TIBDN) from March 2020 to April 2022.Informed consent was obtained and all participating hospitals granted research ethics approval by their respective ethics boards (REB #2218).Demographic and clinical data were obtained from electronic medical record review and patient interview.Serial nasopharyngeal or mid-turbinate samples were collected starting from date of study enrollment until recovery, refusal, or death.
Patient demographics and clinical characteristics were compared between the RDV-treated and untreated groups using Welch's t-test for continuous variables and Pearson's Chi-Square or Fisher's exact tests for categorical variables.Statistical testing and plotting were performed using RStudio v1.4.1717 with the dplyr 27 v1.0.10 and ggplot2 28 v3.3.6 packages.The Mann Whitney U test was used to compare the distribution of mutation counts, and transition and transversion mutation types between patients with and without RDV treatment.The Chi-Square test was used to assess the difference between the Ti/Tv ratios in the RDV-treated and untreated groups.Statistical significance was determined based on a p-value of less than 0.05.

T A B L E 2
Abbreviation: RDV, remdesivir.
demographics and clinical characteristics of cohort.