Evolution of viral variants in remdesivir‐treated and untreated SARS‐CoV‐2‐infected pediatrics patients

Abstract Detailed information on intrahost viral evolution in SARS‐CoV‐2 with and without treatment is limited. Sequential viral loads and deep sequencing of SARS‐CoV‐2 from the upper respiratory tract of nine hospitalized children, three of whom were treated with remdesivir, revealed that remdesivir treatment suppressed viral load in one patient but not in a second infected with an identical strain without any evidence of drug resistance found. Reduced levels of subgenomic RNA during treatment of the second patient, suggest an additional effect of remdesivir on viral replication. Haplotype reconstruction uncovered persistent SARS‐CoV‐2 variant genotypes in four patients. These likely arose from within‐host evolution, although superinfection cannot be excluded in one case. Although our dataset is small, observed sample‐to‐sample heterogeneity in variant frequencies across four of nine patients suggests the presence of discrete viral populations in the lung with incomplete population sampling in diagnostic swabs. Such compartmentalization could compromise the penetration of remdesivir into the lung, limiting the drugs in vivo efficacy, as has been observed in other lung infections.

Here we report the application of similar methods in a personalized medicine approach to investigating the impact of remdesivir on SARS-CoV-2 within an individual and to evaluate potential biomarkers that can be used to monitor clinical efficacy. The data provide further insights into fundamental questions of SARS-CoV-2 evolution and coinfection.

| Sample collection and viral sequencing
Nasopharyngeal swab samples were collected and tested for SARS-CoV-2. Full-length SARS-CoV-2 genome sequences were obtained from all positive samples using SureSelect XT target enrichment and Illumina sequencing. For each patient, a unique patient reference was generated by mapping the remaining reads of the first sample to the SARS-CoV-2 reference genome (NC_045512) from GenBank using bwa-mem. 11 Reads from the subsequent samples of the same patient were mapped to this patient reference. Consensus sequences were aligned using MAFFT. 12 Only genomes with more than 80% genome coverage and a mean read depth of 100 or above were included in downstream analysis.

| Phylogenetic analysis
The maximum likelihood tree of the alignment was constructed using RAxML, 13 with the GTR model and 1000 bootstrap replicates. All trees were rooted on the SARS-CoV-2 reference genome NC_045512.

| Analysis and figure generation
Analysis was completed in R 3.6.1 using Rstudio 1.2. In general, data were processed using the tidyverse family of packages (v1.2.1). We employed the fisher.test in R to compare the count-data of mutations for treated and untreated samples across all individuals. The Mann-Whitney-Wilcoxon test was implemented using the wilcox.test in R. We used Pearson's and Spearman's rank correlation for correlation analysis. This was done using the lm.test and cor.test function in the stats package in R.

| Haplotype reconstruction
Haplotypes were reconstructed using HAplotype Reconstruction Of Longitudinal Deep Sequences (HaROLD) with default settings. 14 HaROLD does not statistically support haplotypes from a single minority variant alleles (MVAs), we therefore constructed by hand the haplotypes for Patient B. In this case, the haplotype frequency was taken to be that of the single MVA.

| Quantification of subgenomic RNA
We employed Periscope to detect subgenomic RNA (sgRNA). 15 sgRNA is identified based on the detection of the leader sequence at the 5ʹend (5ʹ-AACCAACTTTCGATCTCTTGTAGATCTGTTCT-3ʹ) of the sequence. To optimize recovery, we excluded genomes with less than 90% coverage and less than 100 mean read depth (MRD).

| Minority variant calling
Minority allele variants had to have a frequency of above 2% and with a minimum of four supporting reads identified at sites with a read depth of ≥5 using VarScan. 16 Transient MVAs, which occurred at one time point in an individual, were discarded from the analysis.

| Structural biology
The structure of the spike protein PDB 6XR8 was visualized using VMD. Mutations were modeled using the Swiss model.   Table 1. Five patients had pre-existing comorbidities associated with primary or secondary immunodeficiency (Table 1) (Table 1). following the first positive sample ( Figure S2). Of the three patients who received remdesivir, only Patient D had total suppression of viral RNA during treatment followed by a rebound of the virus after treatment cessation (Figure 1). The four ICU patients were clinically most unwell, requiring assisted ventilation. All three remdesivirtreated patients showed clinical improvement after starting the drug, associated with falls in temperature (all) and inflammatory markers (A and G) ( Figure S3). All three were weaned from conventional ventilation before the treatment course was completed with decreases in oxygen requirements. In Patient D, a significant reversal of respiratory deterioration was noted once remdesivir was started;

| Viral load trajectories by cycle threshold and clinical markers of infection
inhaled nitric oxide was stopped within 96 h coincident with weaning from high frequency oscillatory to conventional ventilation. Patient D, who alone required inotropic support, achieved hemodynamic stability off inotropes within 5 days of starting remdesivir.

| Inter-and intrahost phylodynamics of SARS-CoV-2
To investigate the possibility of remdesivir resistance in Patients A, D, and G, we deep sequenced all samples as well as those from untreated patients for comparison ( Figure 1). The sequencing metrics for the samples are summarized in Table S1. Relative to their first available sample, there were nine polymorphisms identified in viruses from Patients A, H, and I, of which five were nonsynonymous with four in the ORF1ab (nsps 1, 3, 4, 5) and one in the Spike protein, S2 subdomain (Table S2). None of the identified SNPs were sites identified as common homoplasies or those known to be susceptible to Illumina sequencing error, none have been associated with remdesivir resistance and none, other than spike mutation P812L, which is located within a predicted CD4 T cell epitope, were in known or predicted immune epitopes. 18-25 P812L was not predicted to alter spike protein structure, although it is not known whether it would abrogate T cell binding ( Figure S4). 18

| No mutagenic signature identified for remdesivir
In vitro studies have not shown lethal mutagenesis to be a feature of remdesivir. 26 [47][48][49] treatment, we compared the mutational burden and patterns of transitions and transversions in treated and untreated patients. In accordance with current understanding of its mode of action, we found neither an increased mutational burden in remdesivir-treated patients nor any evidence of associated mutational signature (- Figures S5 and S6). We found no evidence that the proportion of transitions and transversions distribution of mutation was dependent on the treatment (Fisher t test, p = 0.13).

| Measurement of subgenomic RNA
To determine whether, despite stable viral ct values, remdesivir treatment may have inhibited viral replication we analyzed subgenomic RNA (sgRNA) detected by the genome sequencing using Periscope. 15 Stacked bar-plots of the frequency of sgRNA reads per 100 000 mapped reads (sgRPHT) for each gene and corresponding ct values for each patient are compared in Figure 3A. Unlike previous reports, no correlation was found between sgRPHT and ct-values ( Figure S7A). 28 We have excluded Patient D from the comparative analysis below as their virus was suppressed below the limit of detection during remdesivir treatment.

| Evidence of mixed infection
We observed changes in the consensus sequences at different timepoints in Patients A, H and I ( Figure 2A). As depicted on the tree, no SNPS relative to reference sequence NC_045512 are shared across all patients and neither are there any such SNPs that are unique to patients admitted to ITU and/or on treatment (Table S3) possibly that the mutations they carry are deleterious as has been postulated for influenza. 33 The close clustering between each of the noncirculating viruses with other haplotypes in their cognate hosts supports within-host evolution ( Figure S9), as opposed to coinfection with multiple viruses as has previously been suggested. 34 A possible exception to this conclusion was Patient A. It is conceivable that haplotype Hap1_A, which is identical to the viral sequence from patient D, could have resulted from superinfection.
Patients A and D were colocated in the Intensive Care Unit (ICU) ( Figure S1). Haplotype Hap1_A was not identified in the sample taken from Patient A before their transfer to the ICU but was detected approximately 28 h later in a sample taken 2 h before starting remdesivir ( Figure 4B). However, no other healthcare-associated transmissions were reported in the ICU and local epidemiological investigation suggests it to have been unlikely. Moreover, in view of the rapid appearance of Hap1_A following transfer to ICU and the variation in haplotype frequencies over time, it is more probable that Hap1_A and 2_A were present in the patient's first sample but not detected. This would support the within-host evolution of these three closely related strains, for one of which, Hap3_A, there is no evidence of circulation in the population. The possibility that sgRNA, may, together with viral load, provide a biomarker of response to remdesivir should now be explored.
We observe considerable variation in the consensus SARS-CoV-2 sequences obtained from four of the nine patients, which appears unrelated to remdesivir treatment. The heterogeneity is explained by the presence of between 2 and 4 stable viral haplotypes that vary in abundance ( Figure 4A). From this, we see that, while remdesivir suppressed viral replication in Patient D, the same dose had little impact on the viral RNA from the identical strain, Hap1_A, in Patient A ( Figure 4B).
Remdesivir resistance was not found to account for the difference in  (Figures 1 and 3A). At the same time suppression of viral RNA in Patient D may reflect higher remdesivir tissue levels seen in neonates in which drug clearance is reduced. 41 Taken together with pharmacokinetic-pharmacodynamic (PKPD) models of SARS-CoV-2 viral dynamics, which predict that greater than 90% inhibition of replication is required to interrupt viral replication, the results presented here corroborate predictions that remdesivir (87% inhibition) monotherapy is unlikely to succeed. 42 Based on previous experience of treating serious RNA infections in the lung, we speculate that early treatment with combination therapy may be required for antiviral effect.
The haplotypes we constructed for Patients A, B, H, and I were phylogenetically sister taxa, with generally only one haplotype per patient identical to freely circulating lineages ( Figure S9). The only outlier in this respect is Patient A, which shared a haplotype with Patient D. It is possible that Patient A could have acquired their haplotype through cotransmission or superinfection. 34 43,44 Moreover, evolution of variants in SARS-CoV-2 has also been described in several immunocompromised patients, although in the absence of haplotype reconstruction, whether these represent multiple co-existing genotypes is unclear. 45 We, therefore, conclude that early within-host evolution most logically explains the diversity seen in Patients A, B, H, and I.
A major caveat to our findings is that of our cohort pediatric, for which both the clinical picture and outcome of SARS-CoV-2 infections are known to differ from adults. Notwithstanding, similar patterns of clinical and virological response to remdesivir have been described in adults. 2,7,46 Moreover, both clinical and viral sequence data from the use of repurposed drugs to treat other severe RNA viral infections have shown similarities in adults and children. 10,38 Additionally, although the largest of its kind, our cohort is small, with In summary, we show that treatment with remdesivir is capable of suppressing SARS-CoV-2 viral and subgenomic RNA in vivo and demonstrate that the latter, in particular, needs further investigation as a potential biomarker for monitoring antiviral therapy. Our data suggest that heterogeneous response to remdesivir is not due to resistance but rather is likely to be caused by suboptimal tissue levels.
The patterns of SARS-CoV-2 within-host genetic heterogeneity un-

CONFLICT OF INTERESTS
The authors declare that there are no conflict of interests.