The Fitness of Molnupiravir-Signed SARS-CoV-2 Variants: Imputation Analysis Based on Prescription Counts and Global Initiative on Sharing All Influenza Data Analyses by Country

Molnupiravir was the first oral antiviral approved for SARS-CoV-2 infection at the end of 2021, but its risk-benefit ratio has dramatically decreased in 2022 [1]. On July 18, 2022, a preprint confirmed that molnupiravir increased transition/transversion (Ts/Tv) at days 3–5 in 59 patients recruited for the AGILE CST2 trial (NCT04746183) compared to 65 untreated controls [2]. The results looked comfortable to the authors since there was no treatment-emergent resistance to molnupiravir (no change in the predicted amino acid sequence of NSP12 and NSP14) but that should have been just a part of the risk assessment.

Studies in immunocompetent [3] or immunocompromised [4] mice and hamsters [5] found increased mutations in the SARS-CoV-2 genome. Merck claimed no viable SARS-CoV-2 remained after the 5-day treatment in the phase 2 trials on subjects at standard risk, 89% of whom cleared the virus anyway [6], but no virus isolation was attempted during the phase 3 trial in subjects at higher risk for disease progression (including those immunosuppressed and hence with higher basal viral loads) [7]. In November 2023, a preprint from Merck-funded investigators showed that in 175 out of 194 unvaccinated but immunocompetent participants in the phase 2a RCT, in culture-PCR molnupiravir eliminated infectious virus within 48 h in the nasopharyngeal compartment, left select genome regions completely unaffected (thus identifying regions where mutation likely abrogates infectivity), and did not alter the magnitude or neutralization capacity of the humoral immune response [8].

Concerns about ongoing replication and transmission events in molnupiravir-treated patients who failed to eradicate the virus after molnupiravir were first raised by William A. Haseltine [9, 10], Carl. Bergstrom, and James Hildreth on November 11, 2021, and have been confirmed in a Global Initiative on Sharing All Influenza Data (GISAID) exploratory analysis prepublished in January 2023 [11]. Normally, the Ts/Tv rate in the SARS-CoV-2 genome evolution is about 2:1 [12], while molnupiravir typically induces a 14:1 ratio [13]. In chronically infected patients, >80% of nucleotide changes are nonsynonymous, as a result of strong selective immune pressure on the encoded antigens (with an overrepresentation of mutations affecting the spike antigen [14, 15]).

On October 5, 2022, Ryan Hisner, a science teacher from Indiana active on GitHub as SARS-CoV-2 variant seeker, first reported on Twitter that he had observed in GISAID 4 hypermutated SARS-CoV-2 sequences (BA.5.2.1, BA.2.3, BA.1.1, and BA.4.1.1 [16]) with hallmarks of molnupiravir action, i.e., with a much higher incidence of synonymous transitions homogeneously distributed across the entire genome, mostly from countries with large molnupiravir use (despite after authorization the black market has been widespread in low- and middle-income countries [17]). On November 3, a BM.2 sublineage with a Ts/Tv ratio of 41 was reported from 9 Australian sequences [18], first suggesting that on a large scale apparently unfit mutations can occasionally be fit. Hisner initially proposed 6 criteria to identify molnupiravir-generated variants.

• 1.

  Large number of single-nucleotide variations (SNVs) compared to the consensus sequence

• 2.

  High Ts/Tv ratio, because of the mechanism of action of the drug

• 3.

  Higher-than-usual percentage of synonymous SNVs (not reflecting selective pressures)

• 4.

  Near-random distribution of mutations throughout the genome (again not reflecting selective pressures)

• 5.

  High percentage of C->T and G->A nucleotide mutations

• 6.

  Mutations at “stable” amino acid residues

In this paper, we investigated the GISAID database and exposed a web repository of molnupiravir-signed sequences. We further correlated the number of signed sequences with the number of prescribed molnupiravir courses by country, taking into account sequencing intensity.

We referred each sequence to the Ultrafast Sample placement on Existing tRee (UShER) web portal (https://genome.ucsc.edu/cgi-bin/hgPhyloPlace). Potential molnupiravir-generated SARS-CoV-2 sequences were selected from the entire GISAID database [19] on February 15, 2023, filtering by the following criteria: >13 transitions compared to the closest UShER branch, a percentage of G->A transitions out of the total number of SNVs ≥ 20%, and a Ts/Tv ratio >3. For each GISAID EPI ID, we annotated the closest PANGOLIN-designated phylogeny using Nextclade (https://docs.nextstrain.org/projects/nextclade/en/stable/user/nextclade-cli.html), sampling date, country, number of transitions, number of transversions, number of G->A transitions, number of C->T transitions, number of amino acid mutations, percentage of nonsynonymous SNVs, mean and standard deviation of all SNV positions along the genome.

The numbers of molnupiravir prescriptions for selected countries in a given time interval were systematically searched from web portals of national stringent regulatory authorities (as listed at https://en.wikipedia.org/wiki/List_of_stringent_regulatory_authorities), as well as from PubMed: population, reference period, and number of molnupiravir-signed SARS-CoV-2 sequences per 1,000 prescribed molnupiravir courses were annotated. The sequencing intensity per country in the same timeframe was retrieved from reference [20].

A publicly accessible web portal was generated at https://ukcovid.xyz/molnupiravir, tabulating 1,094 sequences. For each GIDAID EPI ID, a link to UShER tree was provided.

The mutational spectrum mostly consisted of SNVs that were slightly more commonly nonsynonymous than synonymous (60% vs. 40%), excluding immune selection as the main driving force. SNV distribution across the genome was Gaussian across the entire genome, as shown by the average mean and standard deviation of their genome positions (about 14,500 ± 8,500 nucleotides out of a whole genome of 29,000 nucleotides).

Molnupiravir-signed sequences were deposited predominantly from countries where molnupiravir had been authorized (Table 1), although many countries did not report national statistics on molnupiravir prescriptions neither on the specific stringent regulatory authority web portal nor in the peer-reviewed literature. We could find a stable number of molnupiravir-signed sequences deposited in GISAID (0.8 per 1,000 prescribed molnupiravir courses): the highest frequencies (>3 per 1,000) found in England and Scotland can be explained by the >5-fold higher SARS-CoV-2 sequencing intensity in the UK than in other countries (Fig. 2h in ref [20]). Countries that did not authorize molnupiravir at all (e.g., France and Canada) had few to no molnupiravir-signed sequences. We are likely not taking into account patients enrolled in the MOVe-OUT RCT (which, in addition to France and Canada, enrolled at centers in Argentina, Brazil, Chile, Colombia, Egypt, Germany, Guatemala, Italy, Japan, Mexico, Philippines, Russian Federation, South Africa, Spain, Taiwan, UK, Ukraine, and USA [7]).

The mutational spectrum of SARS-CoV-2 is known to have evolved from wild-type to Omicron (e.g., with a 2-fold decrease in the number of G->T transversions) [28], but it cannot explain the molnupiravir signature others [11] and we detected. Our study extends previous observations about the peculiar signature left in the SARS-CoV-2 genome [29], strengthening the imputation of molnupiravir as the driver of our similar viral mutants found in GISAID searches.

The complexity of the molnupiravir-generated quasi-species is likely to be very high, with a plethora of minoritarian sublineages within the swarm that are likely to remain unreported in the majoritarian sequence (which remains the only one finally deposited in GISAID): as such, the real occurrence of molnupiravir-signed sequences is likely an underestimate limited to the very minority of cases where the signature has become dominant at that timepoint. We could retrieve only 1 sequence out of 1,094 (0.45%) with molnupiravir tagged in metadata (EPI_ISL_10910076 from Belgium). Albeit poor data entry could be the main reason for that, it should be considered that SARS-CoV-2 sequences deposited in GISAID are mostly representing patients sampled at diagnosis (before initiation of antiviral treatment), when the viral load is higher and hence the likelihood of succeeding with genome sequencing is higher. Hence, it is logical to assume that molnupiravir-signed sequences in GISAID could represent baseline samples from further generation of transmission rather than from index cases. It makes sense to hypothesize that only the less mutated variants within the quasi-species swarm are still vital and fit enough to be transmitted. These intermediates would exist in accordance with the lethal catastrophe hypothesis [30], but under this scenario final imputation would remain hard to prove, requiring prospective sequencing of strict contacts of molnupiravir-treated subjects. With prescriptions deauthorized or greatly reduced across the globe, such clinical trials would be hard to initiate.

Perno et al. reported that in 8 patients under molnupiravir pressure, SARS-CoV-2 strains were characterized by a 6-fold higher genetic diversity compared to 7 patients under Paxlovid^® pressure and 5 patients under no antiviral pressure, with a peak between day 2 and day 5 [31]. Donovan-Banfield et al. [32] confirmed that the Ts:Tv ratio was significantly increased at days 1, 3, and 5 in 180 subjects within the AGILE-CST-2 RCT (NCT04746183) compared to participants given a placebo. Fountain-Jones et al. [33] reported that 5 immunocompromised patients treated with molnupiravir at 10 days accrued on average 30 new low-mid frequency variants (between 10 and 90% of reads) compared to 4 untreated patients, and that these new mutations could persist and, in some cases, were fixed in the virus population.

Sanderson et al. [11] systematically investigated the GISAID databases for a signature of molnupiravir mutagenesis, finding a specific class of long phylogenetic branches that appear almost exclusively in sequences from 2022, after the introduction of molnupiravir treatment, in countries and age groups with widespread usage of the drug, and in recruits of the AGILE RCT of molnupiravir. On February 23, 2023, EMA recommended against the authorization of molnupiravir, but on March 13 Merck filed an appeal against that decision [34]. Recently, Sanderson released a freely available online tool (http://movbranchapp.streamlit.app) that allows to assess how likely SARS-CoV-2 mutations are to have been selected by molnupiravir using nucleotide contexts. In conclusion, our data suggest that molnupiravir has the potential to generate fit variants that circulate in the general population.

We thank all data contributors, i.e., the authors and their originating laboratories responsible for obtaining the specimens and their submitting laboratories for generating the genetic sequence and metadata and sharing via the GISAID Initiative, on which this research is based.

An ethics statement is not applicable because this study is based exclusively on published literature.

We declare that we have no conflict of interests related to this manuscript.

This research was supported by the Italian Ministry of Health, “Ricerca Corrente - Linea 1 - INMI L. Spallanzani I.R.C.C.S.” funding.

D.M. analyzed the GISAID database, created the public web repository, and revised the manuscript; D.F. retrieved prescription counts by country and wrote the first draft; F.M. revised the manuscript.

The original SARS-CoV-2 sequences are available from GISAID.org. The dataset generated by this study is publicly accessible at https://ukcovid.xyz/molnupiravir/.