Mechanism of reaction of RNA-dependent RNA polymerase from SARS-CoV-2.

We combine molecular dynamics, statistical mechanics, and hybrid quantum mechanics/molecular mechanics simulations to describe mechanistically the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) RNA-dependent RNA polymerase (RdRp). Our study analyzes the binding mode of both natural triphosphate substrates as well as remdesivir triphosphate (the active form of drug), which is bound preferentially over ATP by RdRp while being poorly recognized by human RNA polymerase II (RNA Pol II). A comparison of incorporation rates between natural and antiviral nucleotides shows that remdesivir is incorporated more slowly into the nascent RNA compared with ATP, leading to an RNA duplex that is structurally very similar to an unmodified one, arguing against the hypothesis that remdesivir is a competitive inhibitor of ATP. We characterize the entire mechanism of reaction, finding that viral RdRp is highly processive and displays a higher catalytic rate of incorporation than human RNA Pol II. Overall, our study provides the first detailed explanation of the replication mechanism of RdRp.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the virus that caused coronavirus disease 2019 (COVID-19), emerged in China in December 2019 and rapidly spread over the world, causing a worldwide health threat with more than 5 million fatalities. From a phylogenetic point of view, SARS-CoV-2 belongs to the β genus of the coronavirus family, which includes other highly infective pathogens such as the Middle East respiratory syndrome CoV (MERS-CoV) or the severe acute respiratory syndrome CoV (SARS-CoV). SARS-CoV-2 has a large (30 kb of positive RNA) genome, which forces it to strike a balance between high replication fidelity and genetic diversity.3, 4, 5, 6 A highly processive RNA-dependent RNA polymerase (RdRp) and a proof-reading exonuclease are the crucial elements for maintaining the stability of the viral genome, while at the same time enabling its mutation to adapt to new environments.

RdRp is the core of the replication machinery of the virus and one of the largest proteins in the viral genome (932 residues). It binds to nsp7 and nsp8 to form an active complex,7, 8, 9 one that first uses sense RNA as a template to generate a negative copy and then, in a second cycle, generates new copies of genomic and sub-genomic RNAs. At least two other proteins are involved in the replication process: a 601-residue helicase and a 527-residue proof-reading exonuclease., A simple BLAST query shows that SARS-CoV-2 RdRp is highly conserved within the coronavirus family, but homologs out of this family show a quite low identity (less than 25%), which suggests that we are faced with a quite new protein, but with a surprisingly high functional efficiency.

Due to its central role in the viral infection cycle, RNA polymerases are a major target for fighting RNA viruses.13, 14, 15 Although vaccines began to be available as of November 2021, 50% of world’s population still remains without any dose of it, and the use of effective antivirals is needed to prevent future fast-spreading coronavirus outbreaks due to ineffective immunizations or the emergence of new variants. As of today, the only FDA-approved drug for the treatment of SARS-CoV-2 infection is a RdRp-inhibitor, remdesivir (R), a C-nucleoside (see Figure S1) that was approved for the treatment of Ebola. Being a negative single-strand RNA virus, Ebola is very distant from SARS-CoV-2, but its replication is also dependent on the action of an RdRp.

Several SARS-CoV-2 RdRp cryoelectron microscopy (cryo-EM) structures have been reported since the beginning of the pandemic. The first structure comprised the RdRp complex. Afterward, structures including different RNA duplexes were resolved showing that RNA binding causes no drastic rearrangements in RdRp.,, W. Yin et al. were able to capture a structure where remdesivir has already been incorporated in a nascent RNA strand. Also, a pyrophosphate (PPi) molecule departing from the active site and two Mg2+ near it were found. Recently, structures of RdRp in complex with favipiravir triphosphate (TP) provided new insights of the pre-catalytic state. However, these structures were resolved with the drug in a non-productive conformation or without Mg2+ ions in the active site. So far, however, the limited amount of atomistic data on the mode of binding and reaction mechanism of incorporation of both natural substrates and remdesivir hampers our ability to develop new and more active compounds.

We present here a comprehensive study on the mechanism of action of SARS-CoV-2 RdRp. Atomistic simulations characterize binding interactions and substrate preferences in the active site of the viral polymerase. While remdesivir binds to the viral active site more strongly than its natural counterpart, the opposite is found for human RNA polymerase II. Molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) simulations demonstrate that the enzyme follows a canonical 2-ion reaction mechanism with a catalytic efficiency higher than that of the highly evolved human RNA polymerase II (RNA Pol II). Calculations and biochemical experiments demonstrate that remdesivir triphosphate (RTP; the expected bioactive form of remdesivir) can be recognized and incorporated into nascent RNA with an efficiency only slightly lower than natural nucleotides, i.e., remdesivir is not an inhibitor of nucleotide incorporation. Extended MD simulations failed to detect any dramatic distortion in the RNA duplex due to the presence of remdesivir that would account for its inhibitory properties. Furthermore, no steric clashes were detected when the nascent RNA duplex was displaced along the exit channel (see below), which argues against the hypothesis that steric clashes are responsible for delayed inhibition, suggesting a more specific inhibitory mechanism related perhaps to a transient covalent bond with the enzyme. In summary, our results, obtained through thorough computational simulations, provide the first atomistic description of the mechanism of action for SARS-CoV-2 RNA polymerase and provide clues on the mysterious inhibition mechanism of remdesivir.

Results and discussion

SARS-CoV-2 RdRp active site architecture

The viral protein (see Figure 1A) has a canonical active site, quite similar to those of other polymerases. The two essential Mg2+ are coordinated by the ⍺ and β phosphate groups of the incoming triphosphate nucleotide, as well as Asp618, Asp760, Asp761, Tyr619, and the O3′ terminal of the negative RNA strand (see Figures 1B and 1C), which—according to the circular reaction mechanism for polymerases by de Vivo and coworkers—is expected to be ionized. Classical molecular-interaction-potential calculations, using Mg2+ as the binding cation, further corroborate the two Mg2+ binding sites (see supplemental information and Figure S2). Compared with available structural data,, our catalytically competent active site unveils the slight movement of the Asp618-Tyr619 loop (Figure S2) in order to create well-defined coordination spheres for the two Mg2+ and to be able to proceed with a phosphoryl transfer with high efficiency. Two arginine residues (Arg624 and Arg553) bind the phosphates of the incoming nucleotide, aligning the gamma phosphate for an effective transfer (see Figure 1D). The preferential affinity for ribonucleotide triphosphate substrates (as opposed to deoxyribonucleotide ones) can be explained by the need for north puckering that controls the alignment of the reactive groups, as well as the presence of specific H bonds between the 2′OH group of the nucleoside triphosphate (NTP) and the side chains of Asp623 and Ser683 in the catalytic site (see Figure 1D). Additional hydrogen bonds are found between i+1 2′OH and Ser759 (see Figure 1E). The base specificity is controlled by the complementarity of hydrogen bonding with the template nucleobase (see Figures 1F and S3) and by phosphate coordination, as well as by the residues surrounding the active site that mechanically introduce strong isosteric requirements, altogether making a non-Watson-Crick pairing scheme very unlikely (see Figure 1F). Overall, the structural picture of the active site emerging from EM structures,,, and atomistic simulations strongly suggests that despite a short evolutionary history, SARS-CoV-2 RdRp has all the structural requirements to be an efficient RNA polymerase both in terms of catalysis rate and substrate specificity.

Figure 1

Active site of SARS-CoV-2 RdRp makes it an efficient polymerase

(A) The replication complex of SARS-CoV-2 formed by the nsp12 RdRp enzyme (in green), the nsp8 and nsp7 cofactors (orange and gray, respectively), and an RNA template and nascent strands (yellow).

(B) Scheme depicting the two-metal-ion mechanism used by SARS-CoV-2 RdRp.

(C) The active site contains a well-defined coordination sphere of the two Mg2+ ions.

(D) NTP substrate recognition in the active site of RdRp is mediated by a pair of arginines, an aspartate, and a serine.

(E) The deprotonated 3′ terminal nucleotide is stabilized by one of the catalytic ions in the RdRp’s active site.

(F) RdRp active site pocket enables base specificity between the incoming nucleotide and the template.

Figures were prepared with 3D Protein Imager. A 3D structure representation can be accessed through https://mmb.irbbarcelona.org/3dRS/s/Y17cub.

We also analyzed the SARS-CoV RdRp active site through MD simulations (see experimental procedures and supplemental experimental procedures). SARS-CoV and SARS-CoV-2 RdRp (nsp12) share more than a 96% sequence similarity, and their structures’ root-mean-square deviation (RMSD) is 0.84 Å. In addition, the loops containing the catalytic residues Asp618, Asp760, Asp761, and Tyr619, involved in the coordination of the Mg2+s, are identical in sequence. During our simulations, SARS-CoV RdRp and SARS-CoV-2 RdRp Mg2+ coordination spheres remained almost identical (see Figure S4). The average distances between magnesium ions were found to be 3.9 ± 0.2 and 3.7 ± 0.1 Å for SARS-CoV and SARS-CoV-2, respectively. Also, the distance involved in nucleotide incorporation between Pα an O3′ atoms remained very similar at 3.6 ± 0.2 and 3.5 ± 0.1 Å, respectively. Thus, the observed difference in activity between them should be attributed to other factors.

SARS-CoV-2 RdRp and human RNA Pol II display opposite binding preferences with respect to RTP and ATP

To directly compare the affinity of the binding sites of SARS-CoV-2 RdRp and human RNA Pol II, we performed alchemical free energy simulations, transforming ATP into RTP either in the binding pocket or in free solution. The obtained free-energy differences were directly translated into affinity differences (ΔΔG) using standard thermodynamic cycles as described above. As shown in Figure 2, the viral RdRp is characterized by a sizable preference for RTP, in agreement with available biochemical data, which is a serendipitous fact given that it was originally designed to block an evolutionarily distinct polymerase in the Ebola virus. On the other hand, the lower affinity for the active site of human RNA Pol II suggests that remdesivir will rarely be incorporated into nascent human mRNA, partially explaining the drug’s low toxicity in humans. The observed preferences can be rationalized by inspecting stabilizing interactions inside enzymes’ active sites. The nitrile group of RTP can accept two hydrogen bonds with Thr687 and Asn691 sidechains of RdRp of SARS-CoV-2 (Figure 2B). Similar interactions were observed in previous MD studies., On the contrary, inside RNA Pol II, no hydrogen-bond donors able to stabilize the nitrile group are found (Figure 2C).

Figure 2

Binding preferences in viral and human RNA polymerases

(A) Binding free energies of ATP and RTP inside SARS-CoV-2 RdRp and human RNA Pol II (negative values corresponding to a preference for RTP over ATP).

(B) Most important interactions between RTP nitrile group and RdRp polymerase.

(C) RTP inside human RNA Pol II active site.

The SARS-CoV-2 RdRp activation mechanism

Before investigating the detailed mechanism of nucleotide incorporation, we explored the mechanism of O3′ deprotonation of the 3′-end ribose inside RdRp. This is a mandatory step of nucleophile formation that precedes nucleotide incorporation. We explored whether RdRp is able to activate the 3′-end nucleotide through a self-activated mechanism as has been proposed for other polymerases. As shown in Figure 3 (see also Video S1), once a nucleotide has been incorporated inside RdRp, the newly created PPi molecule is found in perfect arrangement to abstract the O3′ hydrogen atom. Our simulations confirm that one non-bridging oxygen atom of PPi’s γ phosphate group can abstract the proton of the 3′ hydroxyl with a free energy of activation of 8.1 kcal/mol. This barrier is ∼3 kcal/mol higher than the reported value for O3′ deprotonation inside DNA polymerase eta (Pol-η). Interestingly, the arrangement inside RdRp enables the direct transfer of the proton from the O3′ hydroxyl group to the γ phosphate group of the PPi, while in other polymerases, like the aforementioned DNA Pol-η, the β phosphate group firstly deprotonates the O3′ atom and then donates the proton to the γ group. O3′ atom’s deprotonation inside the RdRp step is found to be slightly endothermic as in DNA Pol’s two-metal mechanism, which, as described by us and others, is overcome by subsequent PPi–H release from the polymerase’s active site.29, 30, 31, 32, 33 Finally, following the replication cycle, nucleic-acid translocation and PPi–H departure from RdRp’s active site enable the binding of the subsequent NTP molecule.

Figure 3

Mechanism of activation through O3′ deprotonation inside RdRp

(A) Scheme depicting the activation mechanism used by SARS-CoV-2 RdRp.

(B) Free-energy profile as a function of the proton transfer coordinate (d3-d4), obtained for the proton transfer reaction between the O3′ of the just-incorporated nucleotide and the γ phosphate group of the newly formed PPi.

(C) Transition state (TS) where the H3′ proton belonging to O3′ atom of nucleotide “i” is halfway to being transferred to γ phosphate group of PPi. Reaction coordinate (RC) consisted of the distances between O3′ or Oγ atoms to the H atom to be transferred, displayed as d3 and d4, respectively.

The SARS-CoV-2 RdRp reaction mechanism of nucleotide incorporation

QM/MM simulations were used (see experimental procedures and supplemental experimental procedures) to study the ability of SARS-CoV-2 RdRp to incorporate either a natural triphosphate (NTP; exemplified here by uridine-5′-triphosphate [UTP] and ATP) or RTP (remdesivir-TP) into nascent RNA. In the QM/MM reactant state, the Pα atom of the ATP is 3.6 ± 0.2 Å far from O3′ atom of the terminal nucleotide (see Figure S5). In addition, the O3′ atom is in line to perform the nucleophilic attack on the Pα atom of the ATP, displaying an angle of 173° ± 5° between Oα-Pα-O3′ atoms. While one Mg2+ cation is coordinated by the O3′ atom of the 3′ terminal, the Pα atom of the ATP, and carboxylate atoms of Asp618, Asp760, and Asp761, the other Mg2+ cation is coordinated by the carbonyl group of Tyr619 and carboxylate atoms of Asp618 and Asp760. In addition, the Mg2+-Mg2+ distance is 3.6 ± 0.1 Å. Following the two-metal-ion mechanism scheme, one Mg2+ activates the O3′ atom toward the attack of Pα of NTP, and the other Mg2+ stabilizes the upcoming PPi leaving group, while both are in perfect disposition to stabilize the negatively charged transition state (TS). In the TS, the phosphoryl group is halfway to being transferred to the O3′ atom of the terminal nucleotide (see Figures 4B and S5). The O3′-Pα distance is 1.9 ± 0.1 Å, while the Pα-Oα distance is 2.1 ± 0.1 Å. Also, the Mg2+-Mg2+ distance is slightly reduced by 0.2 ± 0.1 Å, while Mg2+s maintained the same interactions with their respective coordination spheres. Metal-aided nucleotidyl transfer reactions can proceed through associative or dissociative TSs (SN2 or SN1, respectively) depending on the specific enzyme catalyzing the process. It has been observed that while CRISPR-Cas9 proceeds through a concerted associative mechanism, group II introns display a dissociative one. Thus, in light of our results, the reaction inside SARS-CoV-2 RdRp proceeds through a concerted associative TS where the breaking and forming bond lengths are found to be similarly extended. Finally, in the product state, the nucleotide has been fully transfer reflected in a Pα-Oα distance of 3.5 Å, and a PPi molecule is formed (see Figure S5).

Figure 4

RNA elongation inside RdRp of SARS-CoV-2 and human RNA Pol II

(A) Top, free-energy profiles as a function of the phosphoryl transfer coordinate (d1-d2) for the incorporation of a U, an A, or an R to a nascent viral RNA strand. Profile for the incorporation of the inside human RNA Pol II is also shown. The phosphorylation reaction consists of a nucleophilic attack of the O3′ of the terminal nucleotide on the Pα of the triphosphate nucleotide. Bottom, bar plot displaying free energies of activation for the process for different NTPs and enzymes.

(B and C) Active site views of the TSs were found when ATP (B) and RTP (C) (purple C atoms) are the substrates for the elongation reaction inside RdRp. The phosphoryl group is half-way to being transferred. One Mg2+ activates the O3′ toward nucleophilic attack and stabilizes the negatively charged TS. The other Mg2+ stabilizes the charged TS, as well as the nascent negatively charged PPi molecule. Distances involved in the reaction are shown as dotted lines with their average values in Å.

(D) TS insight of the elongation reaction catalyzed by human RNA pol II.

A careful analysis of trajectories and the free-energy profiles shows that the formation of a phosphodiester bond proceeds through a single free-energy maximum (TS) corresponding to an activation barrier of 15.8 kcal mol−1 for UTP and 16.2 kcal mol−1 for ATP (see Figures 4, S5, and S6 and Videos S2 and S3) and is characterized by a negative free energy of ∼5 kcal/mol (prior to PPi release). To obtain an estimate of the relative efficiency of the viral enzyme, we studied the RNA polymerization reaction catalyzed by human RNA Pol II (see Figures 4D and S7 and Video S4). Recently published kinetic rate constants of nucleotide incorporation by SARS-CoV-2’s RdRp provide support for the quantitative accuracy of our estimates, giving confidence to the suggested reaction mechanism. To evaluate the degree of fitness of the enzyme, we used the same procedure to predict the RNA Pol II activation barrier (see Figure 4), getting a value of 17.8 kcal mol−1 for ATP, in good agreement with experimental estimates., This value is ∼1.5 kcal/mol higher than that of the viral RdRp (see Figure 4A), demonstrating that despite its short evolutionary trajectory, the efficiency of SARS-Cov-2 RdRp is at least similar to those of highly evolved eukaryotic polymerases.39, 40, 41

Our equilibrium trajectories demonstrate that RTP fits very well into the active site, showing strong canonical Watson-Crick interactions with the uridine in the template RNA. Being isosteric to adenosine, it achieves a perfect shape complementarity and arrangement of reactive groups in the binding site, predicting that it will act as a substrate rather than an inhibitor. This hypothesis is confirmed by QM/MM simulations showing that incorporation of RTP can happen with a free-energy barrier only slightly larger than that of a natural substrate (10% up to 17.4 kcal mol−1; see Figures 4C and S6 and Video S5). Such an increase is mostly due to a slight misalignment of the Oα-Pα-O3′ attack angle (161° ± 8° for RTP, 173° ± 5° for ATP, and 172° ± 5° for UTP) due to RTP nitrile’s group interactions (see Figures S8 and S9). Thus, we can rule out the possibility of remdesivir inhibiting RdRp by blocking the ATP-binding site of RNA polymerase. On the contrary, our simulation strongly suggests that RTP can be efficiently incorporated in front of uridine in a nascent RNA strand. Again, the order of incorporation predicted by our theoretical calculations, UTP > ATP > RTP, agrees perfectly with recent accurate pre-steady-state kinetic experiments, providing additional support to our calculations (see comparison in Figure S10). Dangerfield et al. also observed that although RTP is incorporated more slowly than ATP inside RdRp, it is incorporated more efficiently than its counterpart due to a higher specificity constant (kcat/Km), which also agrees with our estimates.

Remdesivir is well tolerated in an RNA duplex

Remdesivir is a C-nucleoside with an extra C1′ cyano group, and as such, it might distort the helix, causing delayed inhibition of the enzyme due to hampered displacement of the nascent duplex along the exit channel. To explore this possibility, we performed MD simulations of two RNA duplexes differing only in the substitution of a central r(A·U) pair by r(R·U) one (see experimental procedures). The results, summarized in Figure 5, strongly suggest that remdesivir is well tolerated in an RNA duplex and does not introduce any major structural distortion that would justify termination of RNA synthesis. In particular, there are no significant differences between the hydrogen-bonding stability of r(A·U) and r(R·U) pairs (see Figure 5C), and helical parameters of the duplexes are insensitive to the presence of remdesivir. In summary, our results strongly argue against the idea that a dramatic structural alteration of the RNA duplex is the key factor responsible for remdesivir-induced termination of RNA synthesis. Our results agree with recent cryo-EM structures where remdesivir was incorporated at different or multiple positions of the nascent RNA strand without altering RNA-duplex structure.,

Figure 5

Remdesivir does not distort RNA structure

(A–F) Major (A) and minor (B) groove width, Watson-Crick base pairing (C), roll (D), twist (E), and RMSD (F) of the double helix are not affected when remdesivir is present (purple) with respect to control sequence (black) during MD simulations. The average values across the simulations are shown in black and purple dots for the control and remdesivir-containing sequences, respectively. Average standard deviations are shown as black and purple bars. The RMSD of one of the five replicas is shown in (F) (see also Figure S11).

Remdesivir does not block nascent-strand elongations through steric hindrance

Trying to explore alternative reasons for the inhibitory properties of remdesivir, we slid the nascent RNA duplex along the exit tunnel of RdRp to make the r(R·U) pair by simulating the addition of extra nucleotides, which allowed us to scan interactions of the RNA at several positions along the exit tunnel. After R was incorporated, we were not able to detect any point of steric clash that could justify stopping the polymerase progression (see Figure S12). Interestingly, when three more nucleotides were incorporated, we found Ser861, whose sidechain is located at around 3.7 ± 0.3 Å from the nitrile group of remdesivir (see Figure S13A and S13B). However, considering the flexibility of the Ser sidechain, especially in a well-solvated microenvironment (see Figures S13A, S13C, and S13D), this strongly argues against the hypothesis that nitrile-Ser861 steric clash can explain the inhibitory properties of remdesivir.,, To further discard the steric-clash hypothesis, we performed an alchemical mutational (A→R) scan along the nascent strand, looking at relative preferences of adenosine and remdesivir at positions from (i+2) to (i+5) in a poly-A RNA duplex embedded in the RdRp channel (see Figure 6). Very interestingly, A/R ΔG differences were very similar for positions i+2, i+3, and i+5, suggesting that there are not specific interactions between the incorporated nucleotide and the tunnel residues (as expected for an RNA polymerase designed to have a continuum output flux). On the contrary, significant differences are found at position (i+4), where a local free-energy minimum would hinder the movement of the nascent RNA toward the exit of the tunnel. The stabilization of remdesivir at position i+4 seems related to the formation of direct or water-mediated hydrogen bonds between the cyano-group nitrogen and side-chain hydroxyl of Ser861 (see Figure S13).

Figure 6

Remdesivir elongation along RdRp’s exit channel

(A and B) Alchemically slided remdesivir along the i+2 to i+5 positions of the nascent RNA strand (A) and its associated free energies (B). In the i+4 position, remdesivir is energetically preferred with respect to neighboring positions.

(C) Insight displaying a stabilizing interaction between Ser861 and remdesivir is shown. Alchemically mutated adenines to remdesvir molecules are depicted in cyan, and the two Mg2+ cations of the active site are depicted as pink balls. The Ser861 residue, which is hypothesized to be involved in a steric clash, is shown.

Current biochemical and structural data show that remdesivir is stalled after 3 or 4 nucleotides are incorporated, but the reasons are unclear.,42, 43, 44, 45 Also, biochemical experiments have showed that inhibition can be overcome with higher NTP concentrations., In order to dissect the mechanism, structures of RdRp in post- and pre-translocated states with remdesivir in different positions of the nascent RNA strand have been resolved., A previous MD study also suggested that the delayed termination mechanism could be due to remdesivir’s destabilization of base pair interactions as well as of the nucleotide located in the active site. Other experiments have suggested that when remdesivir is present in the template RNA strand, it would also inhibit RdRp through a secondary mechanism. Some authors have noted as key for main inhibition a steric clash of the nitrile group with Ser861,,43, 44, 45 but our calculations failed to detect any steric clashes, suggesting the opposite, i.e., the vicinities of Ser861 are quite stabilizing for remdesivir. This would suggest that rather than being sterically unable to reach the position, remdesivir-containing RNA might be trapped in this position, hampering further sliding. Analysis of the recognition site (see Figures 6, S13C, and S13D) suggests that Ser861 could even be involved in a water-catalyzed Pinner’s reaction, which would lead to a transient covalent bond of RNA to the enzyme. The lack of stability of the resulting complex47, 48, 49, 50 precludes its experimental detection, but a transient covalent bond could contribute to trap remdesivir at the i+4 position, requiring an extra addition of NTP to escape from this stalling situation.

To further explore the role of Ser861 as a stalling element in the sliding of remdesivir-containing RNA, we explored the exit channels of human RNA polymerases and CoV RdRps (see Figures S14 and S15). In the case of human RNA polymerases, no serine (or similar residue) was found in the expected displacement path of the remdesivir nitrile group (see Figure S14), suggesting that, if incorporated by human polymerases, remdesivir will not stop sliding of the nascent RNA. On the contrary, the same exercise made with other CoV RdRps (see Figure S15) found an amphipathic α-helix with a conserved Ser, which can play the same role of Ser861 in other coronaviruses. We can speculate that the same could happen for Ebola virus, the original targeted virus for remdesivir, but lack of structural information on Ebola RdRp precludes a detailed analysis.

Conclusions

SARS-CoV-2 RdRp is a protein common to other coronaviruses but shows little homology outside the family. As CoVs are just a few thousand years old, the protein has had a limited evolutionary period, and we could expect low efficiency. However, and quite surprisingly, our calculations demonstrate that the enzyme is very efficient, even more than eukaryotic polymerases. The viral enzyme follows a mechanism that is similar to that of bacterial or eukaryotic polymerases with the transferred phosphate being stabilized by 2 Mg2+ ions exquisitely coordinated by acidic residues of the catalytic site, while the phosphates of the incoming nucleotide are stabilized by a network of basic residues. The SARS-Cov-2 RdRp makes use of a self-activated mechanism where the gamma-phosphate group of PPi molecule deprotonates the hydroxylic 3′ terminal, generating the nucleophile that participates in the subsequent incorporation of a nucleotide.

Quite surprisingly, RTP is an excellent substitute of ATP, although it is not well recognized by human RNA Pol II. Simulations also show that after being bound to the viral active site, RTP does not block RdRp, but it is incorporated into the nascent RNA in front of a uridine. The resulting duplex does not show dramatic structural changes, which would hinder displacement of the nascent duplex along the exit channel. In fact, analysis of the displacement of the RNA-containing remdesivir along the exit tunnel of RdRp fails to detect points of steric clashes. Moreover, free-energy calculations show that the i+4 position is quite a favorable site for remdesivir. This suggests that the protein environment around i+4, particularly a hydrated serine, can act as a trap for the nascent RNA duplex, stalling the displacement of the helix by stabilizing remdesivir by either non-covalent or transient-covalent contacts.

Experimental procedures

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Modesto Orozco (modesto.orozco@irbbarcelona.org).

Materials availability

This study did not generate new unique reagents.

RdRp complex setup

The cryo-EM structure of the RdRp (nsp12) of SARS-CoV-2 complexed to nsp8 and nsp7 cofactors and the cryo-EM structure also complexed with an RNA duplex were used as our starting models. The first model was obtained when no other structure was available by aligning the SARS-CoV-2 RdRp with RdRp from a hepatitis C virus X-ray structure, from which RNA, the two catalytic metal ions, and the diphosphate group of an NTP analog were extracted (see supplemental experimental procedures for details). The structure was refined by comparing the cryo-EM structure of RdRp in complex with a full RNA strand (see Figure S16). Finally, we built an ATP, UTP, or RTP molecule inside the active site of RdRp in the expected orientation required for incorporation into the nascent RNA strand. The systems were solvated and neutralized prior to optimization, thermalization (T = 298 K), and equilibration (see supplemental experimental procedures for details). The final equilibrated structures were the starting points for further MD simulations. The SARS-CoV RdRp starting model, including an RNA duplex and an ATP molecule, was built from an available cryo-EM structure and subjected to the same modeling and equilibration protocols as described above. The progression of the polymerization process was simulated by adding additional base pair steps to the RNA duplex, moving the r(R·U) pair from i+1 to i+5 position while keeping constant the reactive alignments at the i position. These structures allowed us to trace the sliding of the nascent RNA duplex along the exit path and check for potential reasons for the R-induced delayed termination of the polymerization reaction.

Human RNA Pol II complex setup

The X-ray structure of human RNA Pol II (PDB: 5FLM) consisting on a polymerase protein complex (composed on 12 subunits), a DNA template, and an RNA transcript were taken and submitted to further modeling. In order to obtain a fully processive protein enzyme, we aligned the aforementioned human RNA Pol II, whose trigger loop is in an opened conformation, to an RNA Pol II X-ray structure (PDB: 2E2H and 2E2J) in the “closed” state. Finally, we extracted one Mg2+ cation and the triphosphate moiety, absents in the X-ray structure (PDB: 5FLM), from the closed-state X-ray structures (PDB: 2E2H and 2E2J) and introduced them to the final model of human RNA Pol II. The final structure was subjected to the same equilibration protocol as described above (see supplemental experimental procedures).

MD simulations on RdRp and human RNA Pol II

Classical trajectories were used to refine and check the stability of complexes prior to running QM/MM simulations, as well as to determine the binding mode and perform free-energy calculations. Minimization, thermalization, and equilibration were performed using standard procedures, as described in supplemental experimental procedures. Production simulations were carried out using the AMBER 19 program and state-of-the-art conditions for a total time of at least 0.5 μs. Water molecules were described through the TIP3P model, parameters of magnesium ions were taken from Allner et al., and Carlson et al. parameters were used for triphosphate groups, PARMBSC1 for DNA, PARMBSC0-chiOL3 for RNA,60, 61, 62, 63 and ff14SB for the proteins. Parameters and charges of ATP, UTP, RTP, remdesivir nucleotide, and 3′ terminal nucleotides and R were derived to be compatible with the force fields making use of the RED server. Additional details of the simulation setups can be found in supplemental experimental procedures and Figure S17.

QM/MM exploration of the minimum free-energy paths and potential of mean force

QM subsystems used in reactivity calculations for RdRp and RNA Pol II are shown in Figure S18. The link-atom method was employed to join QM and MM regions. The hybrid QM/MM models were built using randomly selected snapshots obtained in the last ns of unrestrained MD simulations, which were then minimized and re-equilibrated at a QM/MM-hybrid level of theory.

QM/MM-MD simulations were performed to obtain minimum free-energy paths (MFEPs) by means of the string method., This method allowed us to explore different reaction mechanisms and select the preferred one in terms of free energy. Sixty to one-hundred twenty string nodes were used. Afterward, a path collective variable (CV) was defined to obtain the potential of mean force (PMF) using 60 to 120 umbrella sampling windows. The chosen set of CVs that followed the progress of the reactions and the breaking and forming bonds are shown in Figures 1B and 3A and supplemental experimental procedures (see Figure S18). MFEPs were obtained at the DFTB3,/MM level, and the PMFs were corrected at the B3LYP/6-311++G∗∗ level. Each of the sampling windows consisted on 20 ps of equilibration followed by 200 ps of production. We checked that the length of the production for the PMFs was sufficient to reduce the statistical error to the order of 1 kcal mol−1. The error of all free-energy barriers and profiles was calculated as 95% confidence intervals and reached error values within ±1 kcal mol−1. The AMBER program with electrostatic embedding was used for the QM/MM calculations. Corrections at the high level of theory were made with Gaussian16 (see supplemental experimental procedures).

Binding free-energy calculations

The difference binding free energy of ATP and RTP to polymerases was determined by computing the differences in the free energy associated with the ATP/RTP change in protein-complexed and isolated states. The human and viral polymerases, complexed with an alchemical ATP/RTP residue in the binding site, were simulated in NPT conditions for 250 ns in both physical endpoints corresponding to ATP and RTP, and the initial 50 ns were considered equilibration. As a reference state, a single solvated ATP/RTP residue complexed with Mg2+ was simulated in the same manner. The resulting 200 ns trajectories were then used to extract seeding frames for the non-equilibrium free-energy protocol based on the Crooks theorem. To obtain the binding free energies, 200 short 1 ns simulations were launched in each direction, with λ changing steadily from 0 to 1 or vice versa. Values of non-equilibrium work were computed for each run and converted to free energies using the Bennett acceptance ratio (BAR) method, all using an in-house implementation of the protocol (gitlab.com/KomBioMol/crooks).

Free energy of remdesivir progression

To verify whether steric hindrance could be responsible for the stalling of RNA extension by RdRP, we performed alchemical free-energy simulations in which an adenine was mutated into remdesivir within the poly-A RNA duplex at four different positions (from i+2 to i+5). A non-equilibrium free-energy protocol virtually identical to the one outlined above was used to obtain relative free energies.

MD simulations on nascent RNA

We performed MD simulations on two duplexes, r(CGCGAAUUGCGC)·r(GCGCAAUUCGCG) and r(CGCGARUUGCGC)·r(GCGCAAUUCGCG), to determine the structural impact of the introduction of a remdesivir in a canonical RNA duplex. Starting structures were those expected for a canonical RNA duplex as implemented in AMBER. Systems were hydrated, minimized, thermalized, and equilibrated using standard protocols, before MD were performed at constant temperature (T = 298 K) and pressure (p = 1 atm). A total of 5 μs cumulative simulation time was sampled for each system. Details of simulations are shown in supplemental experimental procedures.