Molnupiravir is an antiviral drug, licensed in some countries for the treatment of coronavirus disease 2019 (COVID-19). In the body, molnupiravir is ultimately converted into a nucleotide analogue, molnupiravir triphosphate (MTP; also known as β-d-N-hydroxycytidine triphosphate)³. MTP is incorporated into RNA during strand synthesis by viral RNA-dependent RNA polymerases, where it can result in errors of sequence fidelity during viral genome replication. These errors result in many viral progeny that are non-viable and thus reduce the virus’s effective rate of growth. Molnupiravir reduced viral replication during 24 h by 880-fold in vitro, and reduced viral load both in animal models¹ and in patients sampled on the final day of treatment². Molnupiravir initially showed some limited efficacy as a treatment for COVID-19 (refs. ^3,4); subsequently, the much larger PANORAMIC trial found that treatment did not reduce hospitalization or death rates in a group of largely vaccinated high-risk individuals². As one of the first orally bioavailable antivirals on the market, molnupiravir was widely adopted in many countries. However, recent trial results and the approval of more efficacious antivirals have since led to several countries recommending against its use^5–7, while longstanding concerns have been raised about potential mutagenic activity in host cells⁸.

MTP is incorporated into nascent RNA primarily by acting as an analogue of cytosine (C), pairing opposite guanine (G) bases (Fig. 1a). However, once incorporated, the molnupiravir (M) base can transition into an alternative tautomeric form that resembles uracil (U) instead. This means that in the next round of strand synthesis, to give the positive-sense SARS-CoV-2 genome, the tautomeric M base pairs with adenine (A), resulting in a G-to-A mutation (Fig. 1b). These G-to-A mutations arise from incorporation of molnupiravir into the negative-sense genome. Incorporation of MTP can also occur during the synthesis of the positive-sense genome: in this scenario, an initial positive-sense C correctly pairs with a G during negative-sense synthesis, but this G then pairs with an M base during positive-sense synthesis. In the next round of replication, this M can then pair with A, which results in a U in the final positive-sense genome, with the overall process producing a C-to-U mutation (Fig. 1c). The free nucleotide MTP is less prone to tautomerization to the oxime form than when incorporated into RNA, and so this directionality of mutations is the most likely⁹. However, it is also possible for some MTP to bind, in place of U, to A bases and undergo the above processes in reverse, causing A-to-G and U-to-C mutations (Fig. 1c).Fig. 1

a, Molnupiravir triphosphate can assume multiple tautomeric forms that resemble different nucleotides. The N-hydroxylamine form resembles cytosine (C); like cytosine, it can pair with guanine (G) while the oxime form more closely resembles uracil (U) and thus can pair with adenine (A) (figure adapted in part from Malone and Campbell⁴⁵). b, In the most common scenario, molnupiravir (M) is incorporated in the N-hydroxylamine form opposite a G nucleotide. It can then tautomerize into the oxime form, which can then pair to an A in subsequent replication, creating a G-to-A mutation. c, Molnupiravir can result in four different mutation types. In the first column, a G-to-A mutation is created by M incorporation opposite a positive-sense G, which can then pair with an A in the next positive-sense synthesis. In the second column, the positive-sense genome contains a C, which results in a G in the negative-sense genome. This G can then undergo the same G-to-A mutation, creating a negative-sense A that finally results in a U in the positive-sense genome, meaning that the entire process results in a C-to-U mutation. Although the biases of tautomeric forms for the free and incorporated MTP nucleotides favour this directionality of mutations, with M incorporated in the N-hydroxylamine form and then transitioning to the oxime form, the reverse can also occur: this results in A-to-G and U-to-C mutations.

It has been suggested that many major SARS-CoV-2 variants emerged from long-term chronic infections. This model explains several peculiarities of variants, such as a general lack of genetic intermediates, rooting with much older sequences, long phylogenetic branch lengths and the level of convergent evolution with known chronic infections^10–13. During the approval process for molnupiravir, concerns were raised about its potential to increase the rate of evolution of variants of concern¹⁴. In response, it was noted that no infectious virus had been isolated at or beyond day five of molnupiravir treatment, and that mutations recovered following treatment were random with no evidence of selection-based bias¹⁵.

During analysis of divergent SARS-CoV-2 sequences, signs of molnupiravir-driven mutagenesis have been noted¹⁶, including indications of possible transmission. Therefore, we aimed to characterize the mutational profile of molnupiravir and examine the extent to which this signature appeared in global sequencing databases.

To establish the mutational profile induced by molnupiravir, we analysed published longitudinal genomic time series that included both untreated patients and patients treated with molnupiravir^17,18, and compared them against a typical SARS-CoV-2 mutational spectrum¹⁹. In agreement with previous findings, we found that molnupiravir treatment led to a substantial increase (8×, 95% confidence interval = 2.9–16.0) in the rate of mutations and that this increase was highly specific to transition mutations (Fig. 2a), especially to G-to-A and C-to-T mutations (hereon we use ‘T’ rather than ‘U’ as in the sequences). While C-to-T mutations are relatively common overall in SARS-CoV-2 evolution^19–21, G-to-A mutations occur much less frequently; therefore, an elevated G-to-A proportion was especially predictive of molnupiravir treatment (Fig. 2b). We looked for evidence of such a signal in global sequencing databases by analysing a mutation-annotated tree derived from McBroome et al.²², containing more than 15 million SARS-CoV-2 sequences from the Global Initiative on Sharing All Influenza Data (GISAID)²³ and the International Nucleotide Sequence Database Collaboration (INSDC) databases²⁴. For each branch of the tree, we counted the number of each substitution class (for example, A-to-T, A-to-G). Filtering this tree to branches involving at least 20 substitutions and plotting the proportion of substitution types revealed a region of this space with higher G-to-A and almost exclusively transition substitutions, which only contained branches of sequences sampled since 2022 (Fig. 2c), suggesting that some change (either biological or technical) had resulted in a new mutational pressure, with mutation classes consistent with those in patients known to be treated with molnupiravir.Fig. 2

a, Comparison of the relative rate of different classes of mutations in typical BA.1 mutations versus those with molnupiravir treatment (molnupiravir data from Alteri et al.¹⁷ and naive data from Ruis et al.¹⁹, scaled to total mutations in naive individuals from Alteri et al.¹⁷; Methods) confirms an elevated rate of transitions, and particularly C-to-T and G-to-A mutations. b, Differences in the proportion of mutations of different mutation classes in individuals treated with molnupiravir (Alteri et al.¹⁷) versus typical BA.1 mutations (Ruis et al.¹⁹) highlight elevated G-to-A proportion as especially indicative of molnupiravir. These are ratios of proportions, so the apparent reduction in transversions does not require an absolute decrease in the number of transversions but can instead be caused by the increased number of transitions. The box plots depict variation over 1,000 bootstrap resamplings, with boxes showing the 25th, 50th and 75th percentiles, and the whiskers having a length 1.5× the interquartile range. c, A scatter plot where each point is a branch with more than 20 mutations, positioned according to the proportion of these mutations that are G-to-A (x axis) or transitions (y axis), reveals a space with elevated G-to-A and transition rate that occurs only with the roll-out of molnupiravir in 2022. d, A change at the same time point is seen when plotting the number of nodes with more than ten mutations and with G-to-A proportion greater than 25%, C-to-T proportion greater than 20% and transition proportion greater than 90%. e, Plotting the number of high G-to-A nodes identified in 2022 against the number of total genomes for each country revealed considerable variation. f, Countries confirmed to have made molnupiravir available had more high G-to-A nodes than countries that did not. The numbers in brackets represent the number of courses of molnupiravir supplied, normalized to population. P = 0.02 for a log-transformed, two-sided t-test. g, Age distribution for US nodes, partitioned according to whether they satisfy the high G-to-A criteria (P < 1 × 10⁻¹⁰, two-sided t-test). Age metadata are missing for some samples, probably non-randomly. Dataset sizes are n = 106 for the high G-to-A branches and n = 2,472 for the other branches. Where a node had many descendants of different ages, age was assigned using a basic heuristic, as described in the Methods. The box plot depicts the minimum, maximum, and the 25th, 50th and 75th percentiles.

We created a criterion for branches of interest, which we refer to as ‘high G-to-A’ branches: we selected branches involving at least ten substitutions, of which more than 25% were G-to-A, more than 20% were C-to-T and more than 90% were transitions. Simulations predicted that this criterion would have a sensitivity of 63% and a specificity of 98.6% for branches involving 13 substitutions (Methods). Branches satisfying the high G-to-A conditions were almost all sampled after the roll-out of molnupiravir in late 2021 and early 2022 (Fig. 2d and Extended Data Fig. 1). Branches were predominantly sampled from a small number of countries, which could not be explained by differences in sequencing efforts (Fig. 2e,f and Extended Data Table 1).

Many countries that exhibited a high proportion of high G-to-A branches used molnupiravir: more than 380,000 prescriptions occurred in Australia by the end of 2022 (ref. ²⁵). More than 30,000 occurred in the UK in the same period^2,26, more than 240,000 in the USA within the early months of 2022 (ref. ²⁷) and more than 600,000 in Japan by October 2022 (ref. ²⁸). Countries with high levels of total sequencing but a low number of G-to-A branches (Canada and France; Fig. 2e,f) have not authorized the prescription of molnupiravir^29,30. Age metadata from the USA showed a significant bias towards samples from older patients for these high G-to-A branches, compared to control branches with similar numbers of mutations but without filtering on substitution type (Fig. 2g). Where age data were available in Australia, they also suggested that high G-to-A branches primarily occurred in an aged population. This is consistent with the prioritized use of molnupiravir to treat older individuals, who are at greater risk from severe infection, in these countries. In Australia, molnupiravir was preplaced in aged care facilities; it was recommended that it should be considered for all residents testing positive for COVID-19 aged 70 years or older, with or without symptoms³¹.

To further probe the link between high G-to-A branches and molnupiravir, we used mutation spectrum analysis, which considers both the types of mutations and the genomic context in which they occur (Extended Data Fig. 4). The spectrum we identified for branches selected by the high G-to-A criteria was, as expected, dominated by transition mutations, with G-to-A and C-to-T directionality more common than A-to-G and T-to-C (Fig. 3a). We also calculated spectra both from patients known to be treated with molnupiravir^17,18 (Fig. 3b) and from general SARS-CoV-2 evolution¹⁹ (Fig. 3c).Fig. 3

a–c, Single-base substitution mutation spectra for high G-to-A branches (a), individuals known to have been treated with molnupiravir (b) and typical BA.1 spectra (c). Each individual bar represents a particular type of mutation in a particular trinucleotide context (Extended Data Fig. 4). Bars are grouped and coloured according to the class of mutation. Within each coloured group, bars are grouped into four groups according to the nucleotide preceding the mutated residue; then, each of these groups contains four bars according to the nucleotide following the mutated residue. The number of mutations has been normalized to the number of times the trinucleotide occurs in the reference genome, and then normalized so that the entire spectrum sums to 1. d, High correlations between spectra from Alteri et al.¹⁷ from patients known to have been treated with molnupiravir, and the spectra from high G-to-A branches identified in this study. Each point represents the normalized proportion of a particular trinucleotide context. Points are coloured so that a context for C-to-T mutations has the same colour as its reverse complement in G-to-A. The values denoted by c are cosine similarity scores.

There was a strong match between the spectrum of known molnupiravir sequences and that of high G-to-A branches, both in terms of the proportions of different mutation classes and of the context preferences within each mutation class (Fig. 3 and Extended Data Fig. 5a). For C-to-T and G-to-A mutations, a comparison of context preferences gave cosine similarities of 0.988 and 0.965, respectively (Fig. 3d). Similar results were seen when comparing to a spectrum calculated from a separate second dataset from a clinical trial of molnupiravir (Extended Data Fig. 5b)¹⁸. The contextual patterns seen in long branches did not correlate with those during typical SARS-CoV-2 evolution (Extended Data Fig. 5c). In high G-to-A branches, G-to-A mutations occurred most commonly in TGT and TGC contexts, which could represent a preference for molnupiravir binding adjacent to particular surrounding nucleotides, a preference of the viral RNA-dependent RNA polymerase to incorporate molnupiravir adjacent to specific nucleotides, or a context-specific effect of the viral proofreading machinery. These correlations between spectra from high G-to-A branches and individuals known to have been treated with molnupiravir strongly support a shared mutational driver, and therefore that the high G-to-A branches are driven by molnupiravir treatment.

Incorporation of molnupiravir during negative strand synthesis results in G-to-A mutations in the virus sequence, while incorporation during positive strand synthesis manifests as C-to-T mutations in the virus sequence after a second round of replication (Fig. 1c). Consistent with this, we observed a strong similarity between the mutational biases in equivalent surrounding contexts within C-to-T and G-to-A mutations when comparing reverse complement triplets with, for example, a G-to-A mutation in the TGC context on one strand being equivalent to a C-to-T mutation in the GCA context on the other strand (cosine similarity = 0.955; Extended Data Fig. 5d).

Although most of the long high G-to-A branches identified had just a single descendant tip sequence in sequencing databases, in some cases we could see that branches had given rise to clusters with a significant number of descendant sequences. For example, a cluster in Australia in August 2022 involved 20 tip sequences, with distinct age metadata confirming they derived from multiple individuals (Fig. 4a). This cluster involved 25 substitutions in the main branch, of which all were transitions, with 44% C-to-T and 36% G-to-A. Closely related outgroups date from July 2022, suggesting that these mutations emerged in a period of 1–2 months. At the typical rate of SARS-CoV-2 evolution, this number of mutations would take years to acquire in an unsampled population with typical dynamics³². Overall in the dataset, we observed a systematic accelerated evolutionary rate in high G-to-A branches (P < 0.001, two-sided t-test), which is consistent with the action of a mutagenic drug.Fig. 4

a, A cluster of 20 individuals emerging from a high G-to-A mutation event. This cluster involves a saltation of 25 mutations occurring within approximately one month, all of which are transition substitutions, with an elevated G-to-A rate. Sequences were annotated with age metadata suggestive of an outbreak in an aged care facility. Phylogenetic placement within the cluster is affected by missing coverage in some regions. b, Examples of further transmission clusters from the UK. Left, four sequences from the UK from February to March 2022 with 13 shared mutations with the high G-to-A signature. Right, a cluster of four sequences from the UK from February 2022 with 31 shared mutations with the high G-to-A signature. c, A sequence from Australia with a high G-to-A signature and a total of 133 mutations relative to the closest outgroup sequence. Just 2 of the 133 mutations observed were transversions; transitions included many G-to-A events. (In the month after this sequence was deposited, two additional related or descendant sequences, EPI_ISL_16315710 and EPI_ISL_16639468, were deposited, which may represent continued sampling from the same patient as they involve a substantial subset of shared mutations, but not full concordance, which is suggestive of complex intrahost evolution).

There are many further examples of high G-to-A branches with multiple descendant sequences, including sequence clusters in the United Kingdom, Japan, USA, New Zealand, Slovakia, Denmark, South Korea and Vietnam (Fig. 4b and Supplementary Information).

During the construction of the daily updated mutation-annotated tree²², samples highly divergent from the existing tree were excluded. This is a necessary step given the technical errors in some SARS-CoV-2 sequencing data, but means that highly divergent molnupiravir-induced sequences might have been missed from this analysis. To search for excluded sequences with a molnupiravir-like pattern of mutations, we processed a full sequence dataset with Nextclade and calculated the proportion of each of the mutation classes among the private mutations (Methods) each sequence carried. This analysis allowed the identification of further mutational events, including some involving up to 130 substitutions (Fig. 4c and Extended Data Fig. 2), with the same signature of elevated G-to-A mutation rates and almost exclusively transition substitutions. The cases we identified with these very high numbers of mutations predominantly involved single sequences and could represent sequences resulting from chronically infected individuals treated with multiple courses of molnupiravir. We verified that nucleotide contexts of the transition mutations observed within the sequence in Fig. 4c were much more likely under the molnupiravir spectrum than the typical BA.1 spectrum (Bayes factor > 10¹⁰).

High G-to-A branches made up a considerable percentage of branches involving more than ten substitutions in some countries (Fig. 5a), suggesting that molnupiravir drives a substantial proportion of large saltations. We found that high G-to-A branches had a different distribution of branch lengths from other types of branches. In typical SARS-CoV-2 evolution, the branch length distribution contains many more nodes with shorter branch lengths than with larger branch lengths; however, for nodes satisfying the high G-to-A criterion this decline was much less pronounced, with long-branch lengths still relatively common (Fig. 5b).Fig. 5

a, Proportion of branches that are high G-to-A for a range of branch lengths in different countries. Data are from collection dates in 2022 and 2023 (submission dates up to June 2023). b, Branch length distributions for high G-to-A and other branches. c, Genomic distribution of mutations in high G-to-A nodes, partitioned into three classes: synonymous mutations; non-synonymous mutations; and non-synonymous mutations that occur four or more times. d, Table of the most recurrent mutations in S in high G-to-A branches. n = the number of high G-to-A branches exhibiting the mutation. Mutation type shows the parental and final nucleotide at the nucleotide position driving the mutation, while context shows the trinucleotide context for the mutated nucleotide, transcribed assuming a NC_045512.2 background.

We also examined whether the mutations identified induced changes to amino acid sequence (non-synonymous mutations) or not (synonymous mutations). We found that for short branches in the tree, 65% (64.6–64.7%) of mutations in the spike (S) gene were non-synonymous. For long branches (ten or more mutations) that lacked the high G-to-A signature, the proportion of S mutations that were non-synonymous was higher at 77% (76.80–77.55%). This increase may reflect, in part, positive selection during intrahost evolution in individuals with chronic infections. In contrast, for long branches with a (molnupiravir-associated) high G-to-A signature, the proportion of mutations in the spike (S) protein that were non-synonymous was 63% (60.3–65.1%), similar to that of short branches, and substantially lower than that of other long branches (P < 0.001).

Despite this overall indication of purifying selection, which is consistent with the actions of a mutagenic drug, there was also evidence for positive selection. Even in high G-to-A branches, there was a concentration of non-synonymous mutations in S, especially among the most recurrent mutations (Fig. 5c). Many of the recurrent S mutations, such as S:P9L, S:A701V, S:K147E, S:R493Q and S:G252S, were also mutations that arise in variants of concern or chronic infections, including S:E340K, which has been associated with sotrovimab resistance (Fig. 5d). There was good correspondence between the contexts in which these mutations occurred and the molnupiravir mutation spectrum. For example, the most common context for G-to-A class mutations among those listed in Fig. 5d is TGT, which has a high enrichment in the molnupiravir spectrum and a low enrichment in the normal BA.1 spectrum (Extended Data Fig. 5c).

There was also a relative concentration of recurrent non-synonymous mutations in the exonuclease encoded by nsp14. This proofreading exonuclease functions to correct errors during genome replication, but typically has poor performance in recognizing mismatches involving molnupiravir³³. Future work could examine whether there is a relationship between specific mutations in nsp14 and tolerance to molnupiravir.

To better test a direct relationship between high G-to-A branches and the use of molnupiravir, linkage analysis was performed for high G-to-A branches sampled in England with treatment data in the Blueteq database³⁴. This analysis found that 31% of clades descended from a high G-to-A branch involved an individual prescribed with molnupiravir (11 sequences in singleton clades, with sampling dates from day 7 to day 61 after the treatment start date). The overall rate of molnupiravir prescription in sequenced individuals in England from 2022 is 0.043%.

Not all branches analysed were linked to a person known to have been prescribed molnupiravir. In some cases, these could represent false positives in our analysis. In addition, the Blueteq database does not contain prescription data for people treated as part of clinical trials (which make up around a third of total molnupiravir prescriptions in the UK) or patients who fell outside interim clinical policy; it is also possible that in some cases we have not sampled the index case treated with molnupiravir, but instead an individual downstream of a treated patient in a transmission chain.

The observation that molnupiravir treatment has left a visible trace in global databases of consensus SARS-CoV-2 genomes, including onward transmission of molnupiravir-derived lineages, will be an important consideration for assessing the effects and evolutionary safety of this drug. Our results are consistent with recent observations in immunocompromised individuals³⁵.

New variants of SARS-CoV-2 are generated through the acquisition of mutations that enhance properties, including immune evasion and intrinsic transmissibility^36,37. The impact of molnupiravir treatment on the trajectory of variant generation and transmission is difficult to predict. A high proportion of induced mutations are likely to be deleterious or neutral, and it is important to consider a counterfactual to molnupiravir treatment that might involve higher viral load, potentially increasing the absolute number of diverse sequences^38,39. However, molnupiravir increases per-sequence diversity in the surviving population, potentially with many mutations per genome, which might provide a broader substrate for selection to act on during intrahost evolution. Importantly, the divergence of the molnupiravir mutation spectrum from standard SARS-CoV-2 mutational dynamics might allow the virus to explore the fitness of distinctive parts of the possible genomic landscape to those that it is already widely exploring in the general population. Molnupiravir-induced mutations could also potentially allow infections to persist for longer by creating a more varied target for the immune system: this might contribute to the PANORAMIC trial’s finding that while individuals treated with molnupiravir had a much lower viral load at day five, they had a slightly higher viral load than placebo-arm individuals at day 14 (ref. ²). Notably, in some countries a significant proportion of sequences with the longest branch lengths are attributable to molnupiravir. However, at the time of writing, the largest clusters satisfying our criteria consist of approximately 20 sequenced individuals.

Considerations of the mechanism of action are important in the design and assessment of antiviral drugs. Molnupiravir’s mode of action is often described using the term ‘error catastrophe’, the concept that there is an upper limit on the mutation rate of a virus beyond which it is unable to maintain self-identity⁴⁰; however, this model has been criticized on its own terms⁴¹ and is particularly problematic in the case of short-term antiviral treatment as it assumes an infinite time horizon. The ‘lethal mutagenesis’ model is much more useful in this context⁴². Not all nucleoside analogue drugs function primarily through mutagenesis. Many act through chain termination^43,44 and therefore would not be expected to cause the effects seen in this study for molnupiravir.

Our study illustrates the far-reaching potential of the extensive genomic dataset created by the community of researchers tracking SARS-CoV-2 evolution. The combination of all available global sequences increased the power of our analyses, while comparisons between countries with different treatment regimens were highly informative. Public health authorities should perform continued investigations into the effects of molnupiravir in viral sequences and the transmissibility of molnupiravir-derived lineages. These data will be useful for ongoing assessments of the risks and benefits of this treatment and may guide the future development of mutagenic agents as antivirals, particularly for viruses with high mutational tolerances, such as coronaviruses.

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-06649-6.