Mutations in the monkeypox virus replication complex: Potential contributing factors to the 2022 outbreak.

Attributes contributing to the current monkeypox virus (MPXV) outbreak remain unknown. It has been established that mutations in viral proteins may alter phenotype and pathogenicity. To assess if mutations in the MPXV DNA replication complex (RC) contribute to the outbreak, we conducted a temporal analysis of available MPXV sequences to identify mutations, generated a DNA replication complex (RC) using structures of related viral and eukaryotic proteins, and structure prediction method AlphaFold. Ten mutations within the RC were identified and mapped onto the RC to infer role of mutations. Two mutations in F8L (RC catalytic subunit), and two in G9R (a processivity factor) were ∼100% prevalent in the 2022 sequences. F8L mutation L108F emerged in 2022, whereas W411L emerged in 2018, and persisted in 2022. L108 is topologically located to enhance DNA binding affinity of F8L. Therefore, mutation L108F can change the fidelity, sensitivity to nucleoside inhibitors, and processivity of F8L. Surface exposed W411L likely affects the binding of regulatory factor(s). G9R mutations S30L and D88 N in G9R emerged in 2022, and may impact the interaction of G9R with E4R (uracil DNA glycosylase). The remaining six mutations that appeared in 2001, reverted to the first (1965 Rotterdam) isolate. Two nucleoside inhibitors brincidofovir and cidofovir have been approved for MPXV treatment. Cidofovir resistance in vaccinia virus is achieved by A314T and A684V mutations. Both A314 and A684 are conserved in MPXV. Therefore, resistance to these drugs in MPXV may arise through similar mechanisms.

Introduction

According to the Centers for Disease Control and Prevention (CDC), the recent Monkeypox virus (MPXV) outbreak has spread to 106 countries (as of September 26, 2022), of which 99 had not previously reported monkeypox (MPX) cases (https://www.cdc.gov/poxvirus/monkeypox/response/2022/world-map.html). A total of 65,415 confirmed cases have been reported as of September 26, 2022 around the world. Due to the rapid spread of MPXV, the FDA has granted emergency use authorization of the Jynneos smallpox vaccine, and drugs tecovirimat, brincidofovir, and cidofovir for MPXV. The Jynneos vaccine, like all smallpox vaccines, is an attenuated vaccinia virus (VACV) [1,2].

MPXV (genus Orthopoxvirus; family Poxviridae) has a double-stranded DNA genome of ∼196,858 base pairs that encodes ∼200 proteins [3], including F8L, which is a family B DNA polymerase (DNA pol). DNA polymerases are critical enzymes for the replication and repair of genomic DNA among various organisms, spanning from archaea to mammals. These enzymes synthesize DNA in template-dependent and -independent manners. DNA polymerases have been classified into seven families: A, B, C, D, X, Y, and RT (reverse transcriptase), based on amino acid conservation and structural homology [[4], [5], [6]]. DNA polymerases that showed sequence homology with E.coli DNA pols I, II, and III were classified into families A, B, and C, respectively. Family D pols are unique to archaea, whereas eukaryotic DNA pols β, λ, μ, and TdT (terminal transferases) belong to family X. DNA pols sharing sequence homology with E.coli pols IV/V belong to family Y. The RT family encompasses DNA pols from retroviruses and eukaryotic telomerase. Regardless of the family, all DNA pols have conserved active site and template-primer (TP) binding modes [7,8]. Family B DNA pols are found in a wide range of organisms (eukaryotes, prokaryotes, bacteriophages, viruses, etc.) [5]. In eukaryotes, at least four DNA pols (α, δ, ε, and ζ) belong to family B. Additionally, viruses such as herpes simplex virus 1 (HSV1), human cytomegalovirus (HCMV), and poxviruses, encode family B DNA pols.

Mutations in DNA pols can alter the phenotype of the virus. For example, a point mutation, N752D, in equine herpes virus 1 (EHV-1) polymerase causes inflammation in the central nervous system and causes poor muscle control leading to gauche movements termed ataxia [9]. In VACV E9 (DNA polymerase catalytic subunit), mutations A498V, A684V, and S851Y in VACV confer an increased mutation frequency in forward-mutagenesis screens [10,11]. The three residues A498, A684 and S951 are conserved between MPXV F8L (DNA polymerase) and VACV E9. To assess if such mutations exist in the MPXV DNA RC, and whether these mutations are contributing to the current outbreak, we conducted a temporal sequence analysis of available MPXV isolates to identify mutations. We leveraged existing data including the crystal structure of VACV E9 protein (family B pol) [12], which shares ∼98% homology with MPXV F8L, to generate a molecular model of MPXV F8L. Additionally, a number of solved structures comprising pols α, δ, ε, ζ, HSV1 pol in complex with TP and inhibitor [[13], [14], [15], [16], [17]], and a state-of-the-art protein structure prediction method, AlphaFold [18], was used to build a DNA RC, and to infer the impact of mutations in 2022 MPXV viruses as well as predict the efficacy of existing vaccines and antivirals in mitigating the 2022 MPXV outbreak.

Materials and methods

Sequence acquisition

The sequences used in the analysis were acquired from the NCBI nucleotide sequence database through the NCBI Virus portal. As of June 28, 2022, all available and complete MPXV nucleotide sequences were downloaded from the portal (n = 205). The sequence dataset included sequences from over 20 countries, during several outbreaks. This included the earliest MPXV sequence in 1965 from Rotterdam, Netherlands, and several sequences from the 2022 outbreak (n = 122).

Sequence processing and analysis

Using Python, an in-house bioinformatics pipeline was developed to rapidly process, and align MPXV sequences [19]. All sequences were compared to a monkeypox virus sequence that was isolated at a Rotterdam Zoo in 1965 (referred to as 1965 isolate hereafter) (NCBI: KJ642614.1). The sequences were cropped to the gene in focus using coordinates provided by NCBI and translated to protein sequences automatically through Expasy [20]. After translation, each sequence was compared with that of the reference protein sequence using the Biopython library [21]. Comparisons were deposited into a data file for further interpretation and the process was applied for the various genes studied (F8L, A22R, E5R, E4R, G9R). The pipeline code is available at GitHub (https://github.com/bluesk1/RapidSequenceAnalysis-MPXV). Using the resultant data files from the pipeline, mutations in 2001 (NC_003310.1), and 2018 (NC_063383.1) were mapped. For the 2022 outbreak, the most prevalent mutations (i.e., present in ≥50% of sequences) in 2022 sequences were included.

Sequences with collection in our dataset were divided into four time periods: 1965–1999 (n = 11), 2000–2009 (n = 26), 2010–2020 (n = 17), and 2021–2022 (n = 126). The sequences from each time period were independently processed through the above-mentioned pipeline to identify mutation prevalence. An in-house circos configuration script (available upon request) was used for the generation of the circos representation of mutations within each time period.

Construction of MPXV replication complex

A low-resolution structure of VACV DNA replication machinery consisting of E9, A20, D4, and D5 has been reported [22]. A PCNA ortholog in VACV has also been identified [23]. Since PCNA is an integral part of the eukaryotic DNA replication machinery, the PCNA ortholog in MPXV (G9R) is likely a part of the MPXV replication fork. Therefore, G9R was added as a fifth component of the MPXV replication machinery, and a minimal MPXV replication fork was constructed in multiple steps as detailed below.

Modeling of the structures of replication complex components

Due to high homology between VACV E9 and MPXV F8L, a homology-derived structure of F8L was constructed by Modeller software [24] using the crystal structure of VACV E9 DNA pol (PDB entry 5N2E) [12] as a template. Similarly, the structure of MPXV E4R in complex with A22R (1-50) was modeled using the crystal structure of the VACV D4/A20 N-terminal (1-50) complex [25]. The structure of the A22R homolog in VACV (A20) is only partially solved. Therefore, we used AlphaFold [18], a state-of-the-art deep-learning molecular modeling program, to generate a molecular model of A22R. To generate a trimeric structure of G9R, we used ColabFold [26], which is based on AlphaFold. Since a high-resolution structure of VACV helicase is not known, and AlphaFold is computationally expensive for generation of a hexamer of VACV/MPXV helicase, we omitted MPXV E5R in our analyses.

Assembly of the RC

To assemble a reliable RC, we first superposed the VACV E9-insert 3 peptide fused to the C-terminal domain of VACV A20 (PDB entry 6ZXP) [27] onto the helix of insert 3 in the modeled structure of MPXV F8L. Next, AlphaFold-generated A22R was superposed onto the A20/E9-insert 3 peptide (PDB entry 6ZXP). Following this step, the crystal structure of the VACV D4 in complex with the A20 N-terminus (PDB entry 4OD8) [25] was superposed on A22R. To this structure, the molecular model of MPXV E4R in complex with the N-terminus of A22R was superposed. Thus, we obtained a complex consisting of F8L, A22R and E4R. To obtain the position of TP and G9R, we superposed the palm subdomain of the S.cerevisiae pol δ holoenzyme consisting of TP, accessory subunits (Pol31 and Pol32), and the PCNA clamp (PDB entry 7KC0) [13] onto the palm subdomain of modeled MPXV F8L. We then superposed the modeled structure of G9R on the PCNA structure. This superposition provided an approximate position of G9R in the MPXV RC. Minor adjustments were made to accommodate for bad contacts of G9R with A22R. A similar guide for the helicase position in VACV, and therefore in MPXV was not available. Hence, MPXV E5R was not included in the MPXV RC.

Results

Mutations in the MPXV RC

We conducted a temporal analysis of available MPXV sequences (n = 204) and identified mutations within MPXV RC components F8L, E5R, A22R, E4R, and G9R (see below) with respect to the 1965 isolate (Table 1 ), which was considered here as a reference isolate. A total of 10 residue positions, where mutations emerged at different time periods were identified (Table 1). Four mutations (2 in F8L, and 2 in G9R) were most prevalent in the 2022 isolates. F8L mutation L108F emerged in 2022, whereas mutation W411L emerged in 2018 and persisted in 2022 isolates. G9R mutations S30L and D88 N emerged in 2022 isolates. Six mutations that appeared in the 2001 outbreak reverted to the 1965 sequence, and revertant mutations persisted in 2022 isolates (Table 1, Fig. 1 ).

Table 1

Mutations in MPXV over various periods.

Gene	Residue	1965	2001	2018	2022	Epitopea
F8LMPXVDNA Polymerase	108	L	L	L	F	ISPDGCYSL (epitope 1)
	411	W	W	L	L	LTFDYVVTF
	428	I	T	I	I	(epitope 2)
	484	A	S	A	A
	501	I	V	I	I
	785	N	D	N	N
A22RProcessivity factor	313	I	S	I	I	–
E5RHelicase (NTPase)	454	D	N	D	D	–
E4RUracil Glycosylase	–	–	–	–	–	–
G9RPCNA ortholog	30	S	S	S	L	-
	88	D	D	D	N	-

Based on Song et al., 2013.

Fig. 1

Circos diagram showing temporal analyses of mutations in the MPXV replication complex: This figure was generated using an in-house Circos configuration script (available upon request) in conjunction with our bioinformatics pipeline in Python (code available at https://github.com/bluesk1/RapidSequenceAnalysis-MPXV). The thickness of the band represents the number of sequences. Each black dot represents a mutation in the MPXV gene.

MPXV RC

The crystal structure of VACV E9 has revealed the topological positions of poxvirus-specific inserts, inserts 0 to 4 [12]. These inserts are also present in MPXV F8L (Fig. 2 a). Compared to the VACV E9 sequence, inserts 0, 1, and 3 are 100% identical, insert 2 has four substitutions, and insert 4 has one substitution. In VACV, only insert 3 has been demonstrated to interact with processivity factor A20 [12]. Therefore, insert 3 in MPXV is also expected to interact with A22R (VACV A20 homolog). Insert 2 is proximal to the ‘fingers’ subdomain (Fig. 2b), suggesting this insert may interact with the fingers as does ‘polymerase associated domain’ (PAD) in Y-family polymerases upon substrate binding [28]. The interactions of inserts 0, 1 and 4 could not be inferred from the modeled RC. The modeled MPXV RC shows that G9R interacts with A22R and the thumb subdomain of F8L. Similar interactions of PCNA (MPXV G9R ortholog) with the thumb subdomain of pol δ have been reported [13].

Fig. 2

Structural details of MPXV DNA pol (F8L) and proposed components of the MPXV replication complex:Panel a shows the ribbon representation of the homology-derived molecular model of F8L. This model was constructed by Modeller [24] using the crystal structure of VACV DNA pol (E9) (PDB entry 5N2E) [12] as a template. Classical polymerase subdomains, namely the palm, fingers, and thumb are colored in red, blue, and green, respectively. The N-terminal domain is colored gray, whereas the 3′-5′ exonuclease (marked as 3′-5′ Exo) is colored olive-green. Poxvirus-specific structures designated as inserts 0, 1, 2, and 3 are shown in (pink, purple, magenta, and deep taupe, respectively). Insert 4 (red) as it is part of the palm subdomain. The template (marked as Temp) and primer (marked as Pri) were derived from the crystal structure of HSV1 (PDB entry 7LUF) [17] by superposing the palm subdomains of HSV1 pol and F8L are colored as cyan (template) and orange (primer). The epitopes of the antibodies induced by VACV-derived vaccines within F8L are marked as Epitope 1 and Epitope 2 (for sequences of these epitopes see Table 1). Epitope 1 (pastel purple) is within the N-terminal domain, whereas Epitope 2 (plum purple) is within the 3′-5′ exonuclease domain. Panel b shows the positions of G9R as derived from the superposition of pol δ in complex with subunits Pol31, Pol32, and PCNA with respect to the structure of F8L. This model was generated by superposing the palm subdomain of pol δ (PDB entry 7KC0) [13] onto the modeled structure of F8L. The PCNA ortholog in MPXV, G9R was constructed using ColabFold [26], a variant of AlphaFold. A22R was constructed using AlphaFold [18] protein structure prediction. The superposition of the G9R molecular model onto PCNA in the pol δ holoenzyme (PDB entry 7KC0), positions G9R close to the thumb subdomain of F8L.

Insert 2 is one of the least conserved sequences between VACV E9 and MPXV F8L with four substitutions within this region. However, insert 2 contains mutations C356Y, G372D and G380S, which confer resistance to phosphonoacetic acid (PAA), a broad-spectrum family B DNA pol inhibitor. C356, G372 and G380 are conserved between VACV E9 and MPXV F8L. These residues are expected to become proximal to the ‘fingers subdomain’ upon ternary complex formation, which induces conformational changes when a polymerase forms a ternary complex consisting of TP, dNTP and divalent cations [[29], [30], [31]]. Thus, PAA resistance mutations may impact the interaction of the ‘fingers’ subdomain with insert 2, and thereby reduce the binding of PAA.

We sought to identify mutations in plausible MPXV RC constituents and conducted structural analyses to understand whether mutations are located at critical positions to influence the function of the replication complex.

Mutation L108F

Mutation L108F was consistently identified in 2022 MPXV isolates. Residue L108 is located in the N-terminal domain of MPXV F8L (or VACV E9) on a loop that is close to the ‘fingers’ subdomain (Fig. 3 a). A rigid-body docking of TP from HSV1 DNA polymerase in complex with DNA (PDB entry 7LUF) and inhibitor PNU-183792 showed that L108 is close to the base moiety of the ‘flipped’ template overhang nucleotide (Fig. 3a). The phenyl group in L108F is closer to the ‘flipped’ base compared to leucine, and would be predicted to increase the binding affinity of polymerase and TP due to increased hydrophobicity. To generate more direct evidence of enhanced TP binding affinity with the L108F mutation, we determined the buried surface area of L108 and F108, and converted this buried surface area into the Gibbs free energy change (ΔG) of binding using the correlation of Chothia [32]. We then used ΔG to calculate equilibrium binding constant in the following equation:

Fig. 3

The location of mutations in MPXV F8L and functional implications:Panel a indicates the position of the L108F mutation in F8L and the interaction of the residue at position 108 with the ‘flipped’ template overhang (rendered as ball-and-sticks). The template nucleotide atoms, carbon, oxygen, nitrogen, and phosphorus, are colored in cyan, red, blue, and orange, respectively. The wild-type residue L108 is shown as brown carbons whereas F108 is shown as gray carbons. This figure also shows the proximity of residue I501, which was mutated to V501 in the 2001 MPXV isolate but reverted to I501 in the 2018 and 2022 isolates. Panel B shows the position of mutations in insert 2 of MPXV F8L. The wild-type residues are colored in magenta carbons, whereas mutations are colored as gray carbons. It is clear from this figure that L411 is exposed to the surface of the protein to have hydrophobic interactions with another protein.

Here, R is the Gas constant, T is temperature in Kelvin, and K is equilibrium dissociation constant. This analysis showed that F108 has ∼2.8-fold greater binding affinity to TP compared to L108. While this difference in binding affinity does not account for the change in entropy compensation, it does suggest that the binding affinity of F108 containing F8L is higher than L108-containing F8L (wild-type). A crucial functional implication of the change in TP binding affinity is that pol processivity is expected to increase [33]. The interactions involving flipped nucleotides with phenylalanine and their biochemical impact have been previously reported in various polymerases [[34], [35], [36], [37], [38], [39]]. It is anticipated that the interaction of F108 with the ‘flipped’ nucleotide will change the biochemical properties of F8L including fidelity, sensitivity to nucleoside inhibitors, and processivity. In different polymerases, these properties were modulated based on the particular mutation and shorter sidechains [37,[39], [40], [41]]. Therefore, mutation L108F in MPXV is a ‘gain-of-function’ mutation, which most likely enhances processivity, changes sensitivity to nucleoside inhibitors, and fidelity of DNA synthesis.

Mutation W411L

W411 resides within F8L insert 2 and is surface exposed, making it ideally positioned to interact with the aforementioned ‘unidentified’ regulatory factor of poxviruses. Hydrophobic residues are typically buried within the protein core; however, F411L would expose a hydrophobic residue to the surface (Fig. 3b). W411L may have evolved to enhance the interaction between insert 2 and a regulatory/processivity factor(s), and the absence of strong electrostatic features of this insert [12] further support our prediction.

Other mutations in F8L

Four mutations in F8L (I428T, A484S, I501V and N785D) emerged in the 2001 outbreak, but reverted to the 1965 isolate sequences in 2018 and 2022. I428 is also in insert 2 and proximal to the PAA resistance mutation (G380S) in VACV E9 [42]. Due to the proximity of insert 2 with the fingers subdomain and the role of the fingers subdomain in nucleotide selection, I428 may be participating in nucleotide selection. Residue S484 is within the 3′-5′ exonuclease domain of F8L and proximal to epitope 1 [43] (Table 1 and Fig. 2). Residue I501 is within the ‘fingers’ domain close to the template nucleotide (Fig. 3a). I501 might function in conjunction with F108 and contribute to the binding of TP to enhance the processivity. Residue N785 is at the bottom of the ‘palm’ subdomain and its functional significance remains unclear.

Mutations in other components of the replication complex

No mutations identified in the E4R (uracil DNA glycosylase) protein. Processivity factor A22R has one mutation that emerged in the 2001 isolate (I313S) but reverted to the 1965 Rotterdam isolate (I313). Similarly, the D454 N mutation in the E5R protein appeared in 2001 and reverted to D454. While a low-resolution structure of VACV (D5), the E5R homolog, has been reported [44], the atomic details and topological position of the D454 equivalent residue in VACV D5 remains unclear. The MPXV G9R homolog in VACV is G8R, which has been reported to be a structural ortholog of human PCNA [23]. There are two mutations in G9R: S30L and D88N that appeared in the 2022 outbreak (Table 1). To understand the function of these residues in MPXV replication, we generated a model of G9R using ColabFold [26]. We superposed the S.cerevisiae pol δ replication complex containing TP [13] and PCNA, which was pre-superposed on the MPXV F8L model (using the palm subdomain). This structure modeling showed that S30 and D88 are in the close vicinity of E4R. Therefore, these mutations may impact the interaction of G9R with E4R. Similar to the F8L mutations L108F and W411L, G9R mutations S30L and D88 N are ∼100% prevalent in 2022 isolates, suggesting a strong correlation among these mutations.

Approved MPXV treatments targeting replication complex

Currently, three antiviral drugs (tecovirimat, brincidofovir, and cidofovir) have been approved for use in MPXV infection. Brincidofovir and cidofovir target the DNA pol (F8L). Brincidofovir is the oral lipid prodrug of cidofovir which must be metabolized into the triphosphate form to be used as a substrate for F8L. Since these drugs are nucleoside analogs, they bind at the dNTP binding site. The triphosphate form of cidofovir [45] has a high affinity with DNA pol and competes with its natural substrate, dCTP. Due to the presence of a 3′OH group, cidofovir does not act as a chain terminator, in contrast to all approved anti-HIV nucleoside inhibitors.

The mechanism of resistance for Brincidofovir and cidofovir has been established. VACV E9 mutations A314T and A684V are cidofovir resistance mutations [45]. Structural modeling and proximity of A314 near the primer of TP suggests that A314T interferes with the shuttling of primer to the exonuclease domain, whereas A684V is in the ‘palm’ subdomain and likely interferes with template positioning, which, in turn, alters the cidofovir-triphosphate binding at the dNTP binding pocket. Therefore, these mutations may have emerged due to cidofovir usage to treat cytomegalovirus-induced retinitis in HIV-infected individuals, who are immunocompromised and resistant to cidofovir.

Conclusions

Of the total 10 mutations in the MPXV RC, 6 reverted to the 1965 isolate. I428 is within insert 2, which also carries PAA resistance mutations but the resistance mutation is not at position 428. It is possible that I428T evolved as a secondary mutation to functionally compensate for the fitness of the viruses that were resistant to PAA (C356T, G372D and G380S). Revertant mutations at positions 484 and 501 are in close proximity to L108. Therefore, it appears that mutations at 484 and 501 reverted to compensate for mutation L108F.

In addition to the mutations in the DNA RC discussed here, other viral proteins from 2022 isolates also contain a variety of mutations. For example, there are four mutations (V6D, D282 N, T601 M, D633E) in the DNA dependent RNA polymerase gene (OPG151) in 2022 isolates. However, due to the mutation in the DNA RC, we have not discussed the mutations in other genes of MPXV. A low-resolution structure of VACV consisting of E9, A29, D4 and D5 has been reported [22]. A PCNA ortholog in VACV (G8R) has also been reported [23]. Reported low-resolution structure does not contain G8R. It is well-known that PCNA and its homologs or orthologs are integral parts of the replication machinery of many species. It is highly likely that VACV G8R and MPXV G9R are also part of the RC in two poxviruses. Our structural modeling shows that G9R can be easily accommodated between F8L and E4R. Additionally, AlphaFold accurately predicted the N-terminal helices of A22R, as seen in the crystal structure of VACV D4 in complex with the A20 N-terminus [25]. These helices were used for the superposition to obtain the position of E4R without any modification in the orientation of the two helices. It is possible that these two helices adopt a different orientation, which may reposition E4R. However, regardless of the position of E4R, there remains enough space in the RC to accommodate G9R within the MPXV G9R.

In summary, we presented the minimum components of the MPXV replication complex and identified mutations in this complex that likely contribute to the 2022 MPX outbreak. Mutation L108F in F8L, which appeared in the 2022 outbreak, is close to the ‘flipped’ template nucleotide and will enhance F8L processivity, whereas W411L likely enhances the binding of a component of the replication complex by altering the surface exposed region. Due to the high sequence and structural similarity among viral B-type DNA pols, the currently-approved antiviral treatments are expected to retain efficacy against the current iteration of MPXV. However, the evolution of resistance mutations remains possible since critical functional pathways have already been susceptible to functional mutations within the viral proteins.

Author contributions

KS conceptualized the study; KS, SRK, and SS wrote the first draft, and KS, SS, SRK, and CLL wrote the final manuscript. SRK, SLK, AR, and KS conducted genetic analyses, and wrote required programs either in R or in Python; KS and SRK conducted structural analysis; KS, SS, SNB, and CLL edited the manuscript and contributed to understanding the pathogenicity of MPXV. All authors approved the final manuscript.

Funding

K. Singh was partially funded by the (Early Concept grant) and the startup support.

Data sharing

Genomic sequences used in this study were obtained from NCBI GenBank.

Institutional review board statement

Not Applicable.

Informed consent statement

Not Applicable.

Declaration of competing interest

CLL is co-founder and Chief Scientific Officer of Shift Pharmaceuticals, Overland Park, KS, USA. KS is a consultant for Sanctum Therapeutics Corporation, Sunnyvale, CA, USA.