Evolutionary and Structural Insights about Potential SARS-CoV-2 Evasion of Nirmatrelvir.

The U.S. FDA approval of PAXLOVID, a combination therapy of nirmatrelvir and ritonavir has significantly boosted our morale in fighting the COVID-19 pandemic. Nirmatrelvir is an inhibitor of the main protease (MPro) of SARS-CoV-2. Since many SARS-CoV-2 variants that resist vaccines and antibodies have emerged, a concern of acquired viral resistance to nirmatrelvir naturally arises. Here, possible mutations in MPro to confer viral evasion of nirmatrelvir are analyzed and discussed from both evolutionary and structural standpoints. The analysis indicates that those mutations will likely reside in the whole aa45-51 helical region and residues including M165, L167, P168, R188, and Q189. Relevant mutations have also been observed in existing SARS-CoV-2 samples. Implications of this analysis to the fight against future drug-resistant viral variants and the development of broad-spectrum antivirals are discussed as well.

Introduction

Coronaviruses (CoVs) are a group of related RNA pathogens.[1] In humans, infections from CoVs can lead to mild to severe diseases.[2] Mild illnesses include some forms of the common cold, while severe varieties comprise severe acute respiratory syndrome (SARS),[3,4] Middle East respiratory syndrome (MERS),[5,6] and coronavirus disease 2019 (COVID-19).[7] Two CoVs were known as human pathogens before 2002. SARS-CoV emerged in 2002. Two more CoVs that are mildly pathogenic were discovered after the SARS pandemic. MERS-CoV, which emerged in 2012, is an addition to this group of human pathogens.[2] Although SARS-CoV and MERS-CoV are highly lethal, their public health, social, and economic damage was relatively small compared to that from SARS-CoV-2,[8] the CoV pathogen of COVID-19. The COVID-19 pandemic has led to catastrophic impacts on the whole world. Its scale and duration are expected to surpass that of the 1918 influenza pandemic.[7,9] As of January 6, 2022, the total number of confirmed global COVID-19 cases was about 300 million, of which 5.4 million succumbed to death (WHO data).[10] Multinational struggles to find therapeutics for both prevention and treatment of COVID-19 have resulted in the authorization of a number of vaccines, antibodies, and small molecule antivirals.[11−14] The global immunization campaigns have produced significant reduction of infections from targeted SARS-CoV-2 variants. However, new SARS-CoV-2 variants that display attenuated response to or complete evasion of vaccine protection have been continuously emerging. The most significant ones were the Delta variant emerging last year and the Omicron variant resulting in the huge infection surge in recently.[15] Although there has not been any solid evidence to show the connection between the use of vaccines and the emergence of new SARS-CoV-2 strains, the widespread use of vaccines has undoubtedly driven significant evolutionary challenges to the virus. For a pandemic at such an unprecedented scale, it is almost inevitable that vaccine-evading strains will evolve due to the selective pressure. So far, most new SARS-CoV-2 strains show mutations in the viral Spike protein, which is directly targeted by vaccines.[16] Luckily, mRNA vaccines allow quick adaption by changing their mRNA sequences to quell newly emerged SARS-CoV-2 strains. However, the race between vaccine development and viral evolution to evade vaccines will likely continue. Since the advent of the pandemic, huge resources and effort have been put into searching currently authorized and investigational drugs for treatment of COVID-19. Although the U.S. Food and Drug Administration (FDA) has approved some repurposed drugs including remdesivir to treat COVID-19, most current evidence shows that these repurposed drugs provide mild benefits to patients.[13] This is disappointing but not necessarily unexpected since many of these repurposed drugs target the human host instead of the virus itself. SARS-CoV-2 has many essential genes that can be directly targeted for antiviral development. These include Spike, RNA-dependent RNA polymerase (RdRp), main protease (MPro, previously called 3C-like protease (3CLPro)), papain-like protease (PLPro), and several other key enzymes.[17] The majority of efforts have been focused on Spike, RdRp, and MPro. Therapeutic antibodies have been developed for Spike. Like vaccines, these antibodies suffer from decreased efficacy toward new SARS-CoV-2 strains.[14] RdRp is the target of remdesivir, which arrests the RNA replication process. Remdesivir is a nucleotide analog.[18] Another nucleotide analog that has been recently approved by the U.S. FDA is molnupiravir.[19] Unlike remdesivir, which inhibits the viral RNA replication, molnupiravir serves as a substrate for the viral RNA replication to introduce extensive mutations that can finally result in nonfunctional virus.[18,19] The mutagenesis consequence of its use needs to be closely monitored. There have been many small molecule inhibitors that have been developed for MPro.[20−38] Nirmatrelvir from Pfizer has been authorized by the U.S. FDA in combined use with ritonavir to treat COVID-19. This combination therapy is branded as PAXLOVID. As a potent inhibitor that directly works on an essential viral enzyme, nirmatrelvir/ritonavir reduced the risk of hospitalization or death by 89% in a phase 2/3 study.[38,39] Nirmatrelvir/ritonavir is being distributed to major hospitals to treat COVID-19 patients and may likely become a game changer for the pandemic. Aside from all the excitement about nirmatrelvir/ritonavir, there remains an ominous question about whether nirmatrelvir/ritonavir will run into the same quagmire that most vaccines are facing right now. In this Perspective, we will try to address this concern from both evolutionary and structural standpoints. YH-53 is another MPro inhibitor that is a promising antiviral to combat SARS-CoV-2 infection.[33] The benzothiazole moiety of YH-53 is involved in direct interactions with multiple active site amino acid residues, aa24–27, that help design potent and specific antivirals. We will discuss it as well. Structures of both nirmatrelvir and YH-53 are shown in Scheme .

Scheme 1

Structures of Nirmatrelvir and YH-53

Conservation of MPro in the Family of Coronaviridae

SARS-CoV-2 is a positive-sense RNA virus.[40] It uses Spike to recognize the human receptor ACE2 for infection.[41] After its entry into a host cell, the genomic RNA is released and translated to two very large polypeptides, pp1a (500 kDa) and pp1ab (800 kDa), in the host cytosol.[42] The processing of pp1a and pp1ab to 16 nonstructural proteins (NSPs) requires proteolytic functions of two internal protease fragments, NSP3 and NSP5, called PLPro and MPro, respectively.[43] Some NSPs assemble into an RNA polymerase complex that replicates both genomic and subgenomic RNAs.[44] Translation of subgenomic RNAs leads to essential structural proteins including Spike for packaging new virions that will be released to infect other host cells. Polypeptides ppla and pp1ab are two separate translation products of ORF1ab, which is the largest open reading frame in the SARS-CoV-2 genome. Their self-cleavage products are absolutely essential for the virus.[45] Of the two internal protease fragments, MPro processes 13 NSPs. Its perturbation is anticipated to have more significant impacts on viral biology than that of PLPro.[17] As such, MPro has been a focus to develop antivirals for SARS-CoV-2.[20−22,46]

SARS-CoV-2 belongs to the sarbecovirus lineage in the Betacoronavirus (β-CoV) genus.[8] Except SARS-CoV and SARS-CoV-2, which are human pathogens, another four viruses in this lineage with published genomes have either bat or pangolin hosts (Table ). Genome sequence alignments of viruses in this lineage indicate 75–90% sequence identity shared with SARS-CoV-2. Nevertheless, MPro proteins from all these viruses share remarkably high amino acid sequence identity, between ∼94% and ∼100%, with that from SARS-CoV-2. MPro proteins from SARS-CoV and SARS-CoV-2 share ∼96% amino acid sequence identity, whereas Spike proteins from the two origins share only ∼76% amino acid sequence identity. Comparing MPro proteins from the sarbecovirus lineage leads to a conclusion that MPro is a highly conserved protein and is far more conserved than Spike. Therefore, MPro-targeting antivirals will be far less of a concern for drug evasion or resistance than Spike-targeting therapeutics. However, this conclusion may not stand well when the comparison is expanded to include other viruses in the β-CoV genus. The β-CoV genus has five lineages including embecovirus, hibecovirus, merbecovirus, nobecovirus, and sarbecovirus.[47] Searching the NCBI database led to 14 published genomes for viruses in the other four lineages. Genome sequence alignments of these viruses over SARS-CoV-2 led to shared identity between ∼52% and ∼57%, which is considerably lower than that shared between sarbecovirus genomes. Extracting MPro proteins from these genomes and aligning their amino acid sequences led to discouraging results about MPro conservation within the β-CoV genus. Except for a hibecovirus that is genetically more related to SARS-CoV-2 than other viruses in these four lineages, MPro proteins from all other viruses display ∼50% sequence identity shared with that from SARS-CoV-2 (Table ).

Table 1

Viruses in the Family Coronaviridae and Their Genome and MPro Amino Acid Sequence Identities Shared with SARS-CoV-2

genus	subgroup lineage	genome accession no.; strain or isolate name	genome sequence identitya shared with SARS-CoV-2 (%)	MPro sequence identitya shared with SARS-CoV-2 (%)
β-CoV	sarbecovirus	NC045512; SARS-CoV-2
		EPI ISL 410721; pangolin CoV P2S	90.31	99.67
		MN996532b; bat CoV RaTG13	96.13	99.35
		EPI ISL 410538; pangolin CoV P4L	85.50	97.06
		NC004718; SARS-CoV	79.70	96.08
		NC14470; bat CoV BM48-31	75.43	94.12
	hibecovirus	NC025217; bat CoV Zhejiang2013	57.24	67.97
	merbecovirus	NC009019; bat CoV HKU4	55.70	50.65
		NC009020; bat CoV HKU5	54.79	51.61
		NC019843; MERS-CoV	55.54	51.61
		NC039207; hedgehog CoV 1	55.93	51.78
	nobecovirus	NC030886; Rousettus bat CoV	54.53	52.61
		NC009021; Rousettus bat CoV HKU9	55.10	52.94
	embecovirus	NC012936; rat CoV Parker	53.13	51.96
		NC001846; murine CoV	52.92	50.98
		NC006213; HCoV-OC43	54.21	48.69
		NC003045; bovine CoV	53.88	48.69
		NC017083; rabbit CoV HKU14	53.81	48.69
		NC026011; β-CoV HKU24	53.35	49.35
		NC006577; HCoV-HKU1	54.42	48.69
α-CoV		NC002306; feline CoV	52.89	44.12
		NC028806; swine CoV	52.55	44.44
		NC038861; gastroenteritis CoV	52.51	44.76
		NC023760; mink CoV	52.90	47.71
		NC030292; ferret CoV	52.39	46.25
		NC032730; Lucheng Rn rat CoV	52.58	41.69
		NC034972; AcCoV-JC34	52.16	42.02
		NC028824; BtRf-α-CoV/YN2012	52.95	45.10
		NC009988; Rhinolophus bat CoV HKU2	52.85	44.44
		NC028811; BtMr-α-CoV/SAX2011	52.82	42.81
		NC028782; camel CoV	46.12	41.04
		NC002645; HCoV-229E	52.48	41.04
		NC005831; HCoV-NL63	53.30	43.97
		NC032107; NL63-related bat CoV	52.34	44.63
		NC022103; bat CoV CDPHE15	52.33	41.69
		NC003436; porcine CoV	52.17	44.95
		NC009657; Scotophilus bat CoV 512	53.02	45.60
		NC028833; BtNv-α-CoV/SC2013	52.47	44.12
		NC010438; Miniopterus bat CoV HKU8	52.12	43.65
		NC010437; bat CoV 1A	52.91	44.30
		NC018871; Rousettus bat CoV HKU10	53.03	44.77
		NC028814; BtRf-α-CoV/HuB2013	52.66	43.79
γ-CoV		NC010646; beluga whale CoV SW1	50.83	42.12
		NC001451; avian CoV	52.54	41.82
		NC010800; turkey CoV	52.61	40.19
δ-CoV		NC016994; night heron CoV HKU19	51.32	32.70
		NC016995; wigeon CoV HKU20	50.85	35.62
		NC016996; common moorhen CoV KHU21	51.41	33.88
		NC039208; porcine CoV HKU15	49.53	35.95
		NC016992; sparrow CoV HKU17	49.59	35.62
		NC011550; munia CoV HKU13	50.61	35.62
		NC016993; magpie-robin CoV HKU18	50.21	35.29
		NC011547; bulbul CoV HKU11	50.91	34.82
		NC011549; thrush CoV HKU12	51.32	35.44
		NC016991; white-eye CoV HKU16	50.26	36.69

The Smith–Waterman algorithm was used to align the genome RNA sequences and MPro amino acid sequences. All alignments were conducted using the SnapGene program.

Bat CoV RaTG13 has three available deposited genome sequences EPI ISL 402131, MN996532.1, and MN996532.2. The most recent deposit MN996532.2 was used in the current analysis. Please note that the EPI ISL 402131 deposit has an MPro amino acid sequence identical to that of the SARS-CoV-2 MPro.

The viral family of Coronaviridae has four genera including Alphacoronavirus (α-CoV), β-CoV, Gammacoronavirus (γ-CoV), and Deltacoronovirus (δ-CoV).[48] α-CoV is more closely related to β-CoV than the other two genera. α-CoVs have 22 complete genome sequence entries in NCBI. Almost all of them share ∼52–53% genome sequence identity with SARS-CoV-2. MPro proteins extracted from these genomes show amino acid sequence identity between ∼41% and ∼47% shared with that from SARS-CoV-2 (Table ). The amino acid sequence identity continues to drop when three γ-CoVs are included in the comparison. Complete genome sequences are available for 10 δ-CoVs. They share ∼50% genome sequence identity with SARS-CoV-2, which is similar to most α-CoVs and γ-CoVs. However, their MPro proteins maintain below 40% sequence identity shared with that from SARS-CoV-2 and the one from night heron CoV HKU9 shares only 32.7% sequence identity with the SARS-CoV-2 MPro. In comparison, Spike proteins from SARS-CoV-2 and night heron CoV HKU9 share 27.9% sequence identity.

Although MPro might be considered generally more conserved than Spike among CoVs, all these sequence alignment results indicate that MPro itself is also a highly mutable protein. The evolution pressure caused by the widespread use of nirmatrelvir/ritonavir and huge number of patients will undoubtedly lead to the emergence of new virus variants with mutations in MPro that cause drug resistance, as previously evidenced by HIV protease inhibitors.[49] Where and what these mutations can potentially be will help guide the development of new antivirals that can be more effective toward new drug-resistant virus variants and prepare us in advance to cope with future viral infection surges. Doing this analysis may also help to design antivirals that are broad-spectrum toward CoVs. We will focus on nirmatrelvir/ritonavir and YH-53 to discuss potential mutations that may emerge in the SARS-CoV-2 MPro.

Crystal Structures of MPro Bound with Nirmatrelvir and YH-53

Nirmatrelvir can be considered as a tripeptidyl drug with a C-terminal nitrile warhead and an N-terminal trifluoroacetamide. It has a P1 β-(S-2-oxopyrrolidin-3-yl)-methyl side chain, a P2 dimethylcyclopropylproline residue and a P3 tert-butylglycine residue (Scheme ). Five crystal structures of the MPro–nirmatrelvir complex (from here on, without specific indication, MPro means the SARS-CoV-2 MPro) have been released to the Protein Data Bank. Two were from the original Pfizer publication about the development of nirmatrelvir that had an early name of PF-07321332.[38] The other three were determined by Zhao et al.,[50] Kneller et al.,[200] and us. Figure A shows an overlay of all five determined structures. As expected, all five structures show close to identical binding modes for nirmatrelvir at the MPro active site. The largest variation is observed at the N-terminal trifluoroacetamide. The oxygen atom of the amide shows the most significant shift of its position in the five structures. Since this amide oxygen is not involved in any hydrogen bonding interactions, this structural shift is likely due to different crystallization conditions that influenced the van der Waals interactions it involves. Proteins in all five structures are largely similar around the active site except at several residues and the small helical region, aa45–51. N119 reorients in all five structures. The nonconserved orientations of N119 have been observed in other MPro–inhibitor complex structures as well. P168, which has van der Waals interactions with the trifluoroacetamide group of nirmatrelvir, also shows large variations in all five structures. Similar observations have also been made in other MPro–inhibitor complex structures. The aa45–51 helix serves as a cap to bind the P2 residue in an MPro substrate. It is a flexible region that has shown weak electron density in some published structures. For all five determined MPro–nirmatrelvir structures, this small helical region shows large variations that lead to very different orientations for residues E47, M49, and L50. In our structure, a water molecule between P168 and nirmatrelvir is clearly observable. A similar water molecule is present in the pdb entry 7VH8 but absent in all other published structures. Nevertheless, the five structures confirm unequivocally how nirmatrelvir interacts with MPro. As shown in Figure B, five protein regions including aa40–44, aa45–51, aa140–146, aa163–169, and aa186–192 interact directly with nirmatrelvir. Mutations in all these regions will likely distort interactions with nirmatrelvir leading to drug resistance. H172 is a critical residue that forms the binding pocket for the P1 residue in a protein substrate. It is viewed as important to maintain the binding to nirmatrelvir. Its mutation will be highly consequential as well.

Figure 1

MPro regions that interact with nirmatrelvir and YH-53. (A) Overlay of five published MPro–nirmatrelvir complex structures. PDB entries are 7RFS (gray), 7RFW (blue), 7SI9 (yellow), 7TE0 (green), and 7VH8 (purple). Residues that display large variations are labeled. (B) Five MPro regions and the residue H172 that interact directly with nirmatrelvir. (C) Six MPro regions and the residue H172 that interact directly with YH-53. The structure is based on PDB entry 7E18.

Two crystal structures of MPro bound with YH-53 have been published.[32,33] They show close to identical binding modes for YH-53 at the MPro active site. As shown in Figure C, YH-53 interacts with the five regions and H172 in MPro that are critical for the binding of nirmatrelvir. In addition, the benzothiazole moiety of YH-53 is also involved in direct interactions with the β strand region, aa24–27. Mutations in these critical regions and residue will likely contribute to viral resistance toward YH-53.

Potential Mutations in MPro to Resist Nirmatrelvir and YH-53 from the Evolutionary Perspective

Mutability of all six regions and H172 in MPro that interact with nirmatrelvir and YH-53 can be partially derived from sequence comparisons between the SARS-CoV-2 MPro and MPro proteins from other CoVs. Figure shows the sequence alignment of MPro proteins from all β-CoVs with all six regions and H172 indicated. Additional alignments of MPro proteins from α-CoVs and β-CoVs and MPro proteins from all CoVs are provided in Figures S1 and S2. Two protein regions aa40–44 and aa140–146 are highly conserved among β-CoVs and among all CoVs. The aa45–51 region displays low conservation among β-CoVs and no conservation among all CoVs. Both aa163–169 and aa186–192 regions are partially conserved among β-CoVs and among all CoVs. H172 is strictly maintained in all CoVs. The aa24–27 region is partially conserved among β-CoVs and among all CoVs. MPro has an N-terminal catalytic domain (NTD) and a C-terminal domain (CTD) that are connected by a short peptide hinge. The peptide hinge displays close to no conservation among all CoVs. In comparison to the NTD, the CTD is less conserved. However, there are residues in the CTD that are strictly maintained among β-CoVs and relatively conserved among all CoVs. Since both the peptide hinge and the CTD do not directly participate in the interactions with nirmatrelvir and YH-53, they will not be discussed further. The C-terminal tetrapeptide in MPro is the recognition site for its self-cleavage. Its conservation among MPro proteins from all CoVs will be informative in assessing the conservation of the MPro active site. As shown in Figure S2, the last residues (P1) of all MPro proteins are strictly maintained as Q among all CoVs, which indicates high conservation around the P1 binding site in the MPro active site. This is corroborated by the sequence alignment results about aa140–146, aa163–169, and H172 that form the P1 binding pocket of MPro and are highly conserved. The second to last residues (P2) vary between four structurally similar residues F/L/M/I indicating some conservation around the P2 binding site in the MPro active site. This partial conservation is supported by the sequence alignment results about aa40–44, aa45–51, and aa186–192. The two regions aa40–44 and aa186–192 are partially conserved, but aa45–51 has no conservation among all CoVs. MPro does not have a P3 binding pocket. This explains that no conservation is observed at the third to last residues among all MPro. MPro has a small P4 binding pocket in the active site that binds its fourth to last residue for self-cleavage. MPro proteins from both SARS-CoV and SARS-CoV-2 prefer a P4 valine in their substrates. Among all CoVs, only V and I are observed as their fourth to last residues with the majority as V indicating relatively high conservation at the P4 binding site. However, the tolerance of I, a residue slightly larger than V, at this site indicates the proneness of this site to mutagenesis as well. Residues that form the P4 binding pocket are from aa163–169 and aa186–189. The partial conservation at these two regions among all CoVs corroborates the partial variability observed at the fourth to last residues of MPro proteins from all CoVs. Therefore, from evolutionary standpoints, MPro mutations in SARS-CoV-2 that are expected to contribute to drug resistance to nirmatrelvir/ritonavir will likely be concentrated in the helical aa45–51 region and distributed in two other regions, aa163–169 and aa186–192, as well. YH-53 involves an additional region, aa24–27, that interacts directly with its benzothiazole. In aa24–27, L27 is strictly conserved among all CoVs. However, all other residues are not conserved. Therefore, in addition to potential mutations in regions aa45–51, aa163–169, and aa186–192, mutations in aa24–27 except at L27 will likely contribute to viral resistance to YH-53.

Figure 2

Sequence alignment of MPro proteins from the Betacoronavirus genus. The Smith–Waterman alignment was conducted using the program SnapGene. The sequence conservation is shown in colored bars where high dark colored bars indicate high conservation, low light colored bars indicate low conservation, and black flat lines indicate close to no conservation. Regions that interact directly with nirmatrelvir and YH-53 are labeled. The hinge that connects the NTD and CTD and the CTD are labeled as well.

Potential Mutations in MPro to Resist Nirmatrelvir and YH-53 from the Structural Perspective

Perceivable mutagenesis in MPro to confer drug resistance can also be partially derived from the structural analysis of MPro bound with nirmatrelvir and YH-53. Figure A shows residues around the P1 binding pocket that bind the nirmatrelvir P1 residue. C145 is the active site cysteine that forms a reversible covalent bond with the nitrile warhead of nirmatrelvir. It is strictly maintained in all CoVs. The nitrile warhead and C145 react to form a covalent thioimidate adduct. The thioimidate nitrogen is within a hydrogen bonding distance to three backbone nitrogen atoms from residues G143, S144, and C145. These three backbone amide amines generate an oxyanion hole that stabilizes the transition state during the hydrolysis of a protein substrate. Due to its preference of an anion bound at this oxyanion hole, the thioimidate nitrogen formed between the nitrile warhead and C145 is likely deprotonated to maintain a negative charge for strong interactions with three backbone amide amines. This may explain the surprisingly high binding affinity of nirmatrelvir toward MPro (Ki = 3.11 nM).[38] In the two pdb entries 7RFW and 7RFS, a proton was added at the thioimidate nitrogen of nirmatrelvir. This may not represent the real situation. Since C145 and the oxyanion hole are critical for the high affinity of nirmatrelvir toward MPro, any change to them will likely distort the binding of nirmatrelvir. However, both C145 and the oxyanion hole are strictly required for the catalysis of MPro. Any subtle change to them may greatly influence the catalysis. The sequence alignment shown in Figure S2 does show that S144 can be mutated to A and T. However, these two are conservative mutations. Therefore, from structural and catalytical standpoints, G143, S144, and C145 are expected to be highly conserved. Their mutations to cause drug resistance are not very likely. Another residue that is critical for the catalysis and activates the catalytic C145 is H41. H41 and C145 form a catalytic dyad.[51] This residue is part of aa40–45 and strictly maintained in all CoVs. Its mutation to generate drug resistance is not expected. N142 is a residue that displays five different orientations in five determined structures of the MPro–nirmatrelvir complex. This residue is partially conserved among all CoVs and can be mutated to C, A, T, or S. Since N142 is part of the binding pocket for the nirmatrelvir P1 residue and is only partially conserved, it is possible that some of its mutations will contribute to viral resistance to nirmatrelvir. Besides N142, several other residues that form the P1 binding pocket include F140, L141, H163, H164, M165, E166, and H172. The lactam amide of nirmatrelvir forms three hydrogen bonds with the F140 backbone oxygen, the H163 imidazole side chain and the E166 side chain carboxylate. The nirmatrelvir P1 side chain mimics Q in a protein substrate. It is expected that the Q side chain amide will form three similar hydrogen bonds as well. These three hydrogen bonds are apparently critical to maintain the enzymatic specificity to recognize Q and digest a protein substrate right after it. Both H163 and E166 are strictly maintained in all CoVs likely due to this reason. Therefore, mutations at H163 and E166 that lead to drug resistance are not likely. F140 uses its backbone to interact with the nirmatrelvir P1 residue and probably the P1 Q in a protein substrate. However, its aromatic side chain forms a π–π stacking interaction with H163 to juxtapose the H163 imidazole side chain to interact with the P1 Q side-chain amide. This interaction is apparently critical to MPro for its substrate selectivity since F140 is strictly maintained in all CoVs. Another residue that is relatively far from the active site but important in juxtaposing H163 to interact with the P1 Q residue in a protein substrate is Y161. Y161 forms a hydrogen bond with the H163 side chain imidazole. It is strictly maintained in all CoVs as well. Based on all the preceding analysis, we can conclude that mutations in F140, Y161, H163, and E166 to cause drug resistance have low probability due to their critical roles in determining the substrate selectivity of MPro. The backbone oxygen of H164 forms a hydrogen bond with and the backbone of M165 is involved in van der Waals interactions with the nirmatrelvir P1 residue. Their mutations will have low influence on the P1 binding site. H163 can be mutated to Q in some CoVs. Q is functionally similar to H. The influence of an H163Q mutation on the binding of nirmatrelvir will not be significant. M165 is a major residue that forms the P2 and P4 binding pockets. We will discuss it later. Another residue that uses its backbone to interact with the nirmatrelvir P1 residue is L141. This residue is partially conserved among CoVs among L, I, V, and M. Its mutations are likely but expected not to be detrimental to the binding of nirmatrelvir. Based on the above analysis, we can predict that mutations around the P1 binding pocket that likely contribute to viral resistance to nirmatrelvir will be at the residue N142. Since this site is so close to the catalytic center, tolerable mutations cannot influence the catalysis dramatically. Therefore, the impact of N142 mutations on the binding of nirmatrelvir is not expected to be dramatic either. Mutations at L141, S144, and H164 are likely but are expected not to be major contributors to drug resistance.

Figure 3

Residues involved in direct interactions between MPro and two ligands, nirmatrelvir and S-217622. (A) MPro residues involved in the recognition of the P1 residue and the covalent warhead of nirmatrelvir. (B) MPro residues that interact with the P2 residue of nirmatrelvir. (C) MPro residues that interact with the P3 residue and the N-trifluoroacetamide cap of nirmatrelvir. Images in panels A–C are based on PDB entry 7TE0. (D) MPro residues that form a recognition pocket for the benzothiazole compound YH-53.

As shown in Figure B, amino acid side chains that contribute to the formation of the P2 binding pocket are those of H41, M49, M165, D187, and Q189. The side chain of M165 interacts also directly with the N-terminal trifluoroacetamide of nirmatrelvir. Among CoVs, M165 varies among M, I, and L. These are similarly functional residues. Although they are not expected to influence the catalysis and substrate selectivity of MPro significantly, mutating M165 to I or L will have dramatic impacts on the binding of nirmatrelvir. Therefore, mutations at M165, most likely to I or L, will lead to viral resistance to nirmatrelvir. They are anticipated to be observed. We have discussed H41. It is part of the catalytic dyad and not expected to undergo mutations. Both the backbone and side chain of D187 form van der Waals interactions with the nirmatrelvir P2 residue. Any change will potentially influence the binding of nirmatrelvir. However, D187 has been strictly maintained in all CoVs. D187 forms a salt bridge with R40, which is another strictly maintained residue in all CoVs. Both residues and their interactions are likely essential to maintain the structural integrity of the enzyme. Therefore, mutations at them that lead to drug resistance are not expected. The side chain of Q189 forms van der Waals interactions with nirmatrelvir at its P2, P3, and N-terminal trifluoroacetamide positions. This residue has low conservation among all CoVs and varies among Q, E, A, and P. Mutations of Q189 to a small residue such as A or P will likely diminish the binding of MPro toward nirmatrelvir. From both structural and evolutionary perspectives, mutations at Q189 to smaller residues that lead to viral resistance to nirmatrelvir are expected. M49 uses its whole side chain to interact with the nirmatrelvir P2 residue and can be viewed as functionally important for the binding of nirmatrelvir. However, it is part of the aa45–51 that shows no conservation among all CoVs. The aa45–51 region shows weak electron density in some of the determined MPro structures indicating its structural flexibility. This may explain its apparent lack of conservation among all CoVs. MPro proteins from different CoVs also show different peptide chain lengths around the aa45–51 region (Figure S2). All these observations mean that the aa45–51 region is highly mutable, and some deletion or insertion may also be likely. Therefore, mutations in the aa45–51 region that contribute to strong viral resistance to nirmatrelvir are expected. Since this is a highly flexible and mutable region prone to deletion and insertion, it is difficult to predict exact mutations that could be expected. R188 uses its backbone to interact with the nirmatrelvir P2 side chain. R188 is not conserved among all CoVs and varies among R, K, A, Q, and E. R188 is a positively charged residue. Its mutations to residues such as E may significantly influence the binding of nirmatrelvir and lead to drug resistance. Based on all analysis related to residues around the P2 binding pocket, we may predict that mutations at M165, R188, Q189, and the whole aa45–51 region will lead to viral resistance to nirmatrelvir.

The P4 binding pocket as shown in Figure C, which binds the nirmatrelvir N-terminal trifluoroacetamide, involves side chains of M165, L167, P168, Q189, and Q192 and backbones of E166 and T190. The potential impacts of M165 mutations have been previously discussed. E166 uses a backbone oxygen to form a hydrogen bond with the trifluoroacetamide nitrogen of nirmatrelvir. Since E166 is strictly conserved, it will not be discussed further. L167 is partially conserved. The only other variation at this site among all CoVs is F, a functionally similar residue. However, an F mutation will likely introduce a large structural rearrangement around the P4 binding pocket that distorts the binding of the nirmatrelvir N-terminal trifluoroacetamide since the side chain of F is significantly bigger than that of L. Therefore, mutations at L167 will likely lead to viral resistance to nirmatrelvir. The P168 side chain forms van der Waals interactions with the nirmatrelvir N-terminal trifluoroacetamide. These interactions are mainly due to the secondary amino acid nature of P168. Among all CoVs, P168 shows no conservation and varies among P, S, A, G, and N. Mutating P168 to a primary amino acid will potentially remove its van der Waals interactions with nirmatrelvir leading to drug resistance. Q189 has been discussed previously and will not be repeated here. Q192 is a strictly maintained residue among all CoVs. It is not expected to undergo mutations to introduce drug resistance. Although T190 can vary among several small residues among all CoVs, these mutations are not expected to dramatically change the binding of nirmatrelvir. So its mutations are likely to be observed but not perceived as major contributions to drug resistance to nirmatrelvir. Therefore, mutations around the P4 binding pocket that likely contribute to viral resistance to nirmatrelvir will be at M165, L167, P168, and Q189.

YH-53 involves similar interactions with the five regions and H172 in MPro that bind nirmatrelvir. Potential mutations in these regions and residue that lead to viral resistance to YH-53 will not be discussed further. YH-53 contains a benzothiazole beyond its C-terminal ketone warhead. As shown in Figure D, the binding of YH-53 induces the formation of a binding pocket for the benzothiazole moiety. This pocket involves residues including H41, M49, T25, and L27. The benzothiazole nitrogen also interacts with the N142 backbone nitrogen through two hydrogen bonds that are bridged by a water molecule. H41 is strictly conserved as we discussed previously. M49 is in the highly mutable aa45–51 region. L27 is also strictly maintained among all CoVs. However, T25 varies among T, M, S, N, V, and I. It has no conservation among all CoVs. The side chain of T25 forms van der Waals interactions with the benzothiazole. Its mutation to S or N will remove these interactions, although impacts of other mutations are hard to predict. Therefore, mutations that likely lead to deteriorated binding to YH-53 around the benzothiazole binding pocket will be potentially observed at T25 and M49.

So based on all analysis, we speculate that mutations in MPro that will potentially lead to viral resistance to nirmatrelvir will likely include residues in the whole aa45–51 region and residues including M165, L167, P168, R188, and Q189. Other mutations that will influence the binding of nirmatrelvir but are not expected as major contributing factors are include those at residues L141, N142, S144, and H164. For YH-53, mutations at T25 in addition to most sites of MPro that were discussed for nirmatrelvir are expected to lead to drug resistance. In the original Pfizer paper, Ki values of nirmatrelvir toward several CoV MPro proteins were reported. In comparison to a Ki value of 3.11 nM for the SARS-CoV-2 MPro, nirmatrelvir has Ki values of 226, 187, and 189 nM toward NL63-CoV (NC005831), MERS-CoV (NC019843), and HKU1-CoV (NC006577), respectively.[38] All three CoV MPro proteins show very low sequence similarity with the SARS-CoV-2 MPro at regions aa45–51 and aa186–192, supporting conclusions drawn above. At the region aa163–169, all three proteins also display 1 or 2 variations from the SARS-CoV-2 MPro. Therefore, the widespread use of nirmatrelvir will likely trigger mutations at these three regions, leading to drug resistance. Nirmatrelvir may be experimentally used to challenge the current SARS-CoV-2 strains to develop drug-resistant strains that can be analyzed to identify mutations in MPro. This may allow matching with the analysis above and facilitate drug design to combat drug-resistant strains. However, stringent safety protocols to eradicate all drug resistant strains have to be strictly followed to avoid spreading them out of research laboratories. Therefore, precautions need to be seriously taken before these experiments are conducted.

Other Concerns

So far, the analysis has been primarily focused on active site residues that directly interact with nirmatrelvir. It is well known that residues in an enzyme that are distant from the active site can influence the active site through coordination with other residues.[52,53] Although it is beyond the scope of this perspective, mutations that are distant to the MPro active site and their combinations with other mutations to cause drug resistance are expected to be observed. Nirmatrelvir contains a reactive warhead toward activated serine and cysteine residues. In the original Pfizer study, the selectivity of nirmatrelvir over several human proteases including caspase 2, chymotrypsin, elastase, thrombin, and cathepsins B, D, and L was tested.[38] It showed that nirmatrelvir has high selectivity over these enzymes. However, the scope of this study is relatively small. Due to the high reactivity of the nitrile warhead, we cannot rule out the possibility that nirmatrelvir inhibits some human proteases as well. It has been shown that multiple human proteases serve critical functions in the SARS-CoV-2 life cycle.[54] Although inhibiting these enzymes provides advantages to contain the virus, this additional antiviral mechanism will also leave the virus with additional routes to achieve drug resistance that might not necessarily involve MPro mutagenesis and consequently complicate the investigation of drug resistance mechanisms. So more stringent characterizations of nirmatrelvir on its selectivity over human serine and cysteine proteases are necessary.

MPro Mutations from Sequenced SARS-CoV-2 Genomes

Global Initiative on Sharing Avian Influenza Data (GISAID) has a large database of sequenced SARS-CoV-2 genomes from collected samples all over the world. From this database, 6 242 921 full MPro sequences without ambiguous sequence assignments can be extracted. Among them, 800 721 mutations can be identified. Although baseline mutations might be due to errors in the sequencing process, significant mutations are observable at many MPro residues as shown in Figure A. Three residues, L89, K90, and P132, have the most mutations. Mutations at these three residues are not expected to produce drug resistance toward nirmatrelvir. Their formations are likely due to virus fitness. Notably, a recent report has found that the P132H mutation in MPro of Omicron strain decreases the thermal stability of MPro.[55] P132 is the first residue of the hinge area. The sequence alignment has demonstrates the it is nonconservative among coronavirus strains. Among regions and sites that interact directly with nirmatrelvir and YH-53 that are shown in Figure B,C, aa40–44 and aa140–146 show low mutation rates, which agrees with conclusions that were drawn based on the sequence alignment among all CoVs. In aa40–44, I43 shows the most mutability, which also matches the sequence alignment of all CoVs at this site. Within the aa140–146 region, N142 shows the most mutability. This agrees well with the sequence alignment results. The top two observed mutations are N142S and N142D. Both mutations are likely to decrease the binding of nirmatrelvir and YH-53. Another region that shows low mutability is aa163–169. But M165, P168, and T169 display much higher mutability than other residues in this region. The top two mutations at M165 are I and L, which is exactly what the sequence alignment results predicted. Both mutations are expected to significantly influence the binding of both nirmatrelvir and YH-53. The top mutation at P168 is S. This mutation is expected to influence significantly the binding to nirmatrelvir as well. As the sequence alignment results predicted, H172 shows very low mutability. Regions that display high mutability are aa24–27, aa45–51, and aa186–192. Both aa45–51 and aa186–192 have higher mutability than aa24–27. The high mutability of aa45–51 was predicted by the sequence alignment results. The most abundant mutation at M49 is I. How exactly this mutation will change the binding of nirmatrelvir and YH-53 needs to be assessed. Since the aa45–51 region is structurally flexible, mutations at other sites, four of which have high occurrence of mutations, are expected to influence the binding of nirmatrelvir and YH-53 as well. Although the whole aa186–192 region shows high mutability, R188 and Q189 have lower mutability than other residues in this region. But mutations at R188 are clearly higher than the baseline for the whole protein. For the aa24–27 region, both T24 and T25 show high mutability. T25, which interacts directly with the benzothiazole moiety of YH-53, has the top two mutations as I and A. Both mutations will greatly influence the binding of YH-53. Since all of these mutations have occurred naturally without evolution challenges imposed by an MPro inhibitor, there is no question that more concentrated mutations will show up in the presence of an inhibitor and potential drug-resistant SARS-CoV-2 variants will develop. It should be noted that those mutants without significant penalty to viral fitness will be more likely to develop drug-resistant variants. Although MPro mutations that have been detected so far are distressing, they can certainly help us prepare for and cope with the emergence of drug-resistant SARS-CoV-2 variants.

Figure 4

Mutation counts at different MPro sites from sequenced SARS-CoV-2 genomes that have been deposited in GISAID; 800 721 mutations were identified in 6 242 921 sequenced MPro entries. (A) Distribution of mutations on all MPro sites. Three sites with most mutations are labeled. (B) Distribution of mutations on aa1–60. (C) Distribution of mutations on aa135–195. In panels B and C, MPro regions and sites that interact directly with nirmatrelvir and YH-53 are colored in red.

Implications

Nirmatrelvir uses its P2 residue and the N-terminal trifluoroacetamide to bind tightly to the P2 and P4 binding pockets of MPro, which are prone to mutagenesis. To better prepare us to contain new virus variants that develop resistance to nirmatrelvir, possible mutations may be introduced to MPro on purpose to study the drug resistance mechanism and use the obtained information to guide the development of effective inhibitors for mutated MPro variants. For example, M165I and M165L might be introduced to MPro to assess their influence on binding to nirmatrelvir and the viral resistance to nirmatrelvir as well. MPro has a highly conserved P1 binding pocket and aa40–44 region among all CoVs. In order to develop antivirals that have potentially lower drug resistance concerns and can potentially be broad-spectrum against all CoVs, one may think to focus on targeting the P1 binding pocket and sites around the aa40–44 region. Although the aa45–51 region is a highly mutable and flexible region, the P2 binding pocket itself is relatively conserved in its recognition of F, L, M, and I. Small molecules that engage only the P1 and P2 binding pockets will likely achieve lower drug resistance and engender broad-spectrum activity against all CoVs. In addition, new technology such as proteolysis targeting chimera (PROTAC) will provide an alternative solution to overcoming drug resistance. PROTACs typically show greater advantages in drug resistance compared to traditional small-molecule inhibitors because (i) they can recruit targets via any binding site where functional and sustained inhibition is not required; (ii) weak binders can become potent degraders; (iii) they degrade the whole target protein; (iv) they have high efficacy due to catalytic mode of action; and (v) they have enhanced specificity due to multiple mechanistic steps.[56−61] Based on the emerging PROTAC technology and existing MPro inhibitors, the development of small-molecule degraders for specific degradation of MPro could be a novel strategy for developing next-generation anticoronavirus drugs with higher resistance barriers.[23,62−64] Although it has many potential advantages, antiviral PROTAC development will face many challenges including pharmacokinetic challenges such as poor cell permeability and low bioavailability associated with high molecular weights and only a few E3 ubiquitin ligases that can be used in PROTAC design. With tremendous effort and fast progress that has been made in the field of PROTAC in drug discovery, it is worth believing that PROTAC technology will provide a promising approach to develop next-generation antivirals against coronaviruses.