Phylogenetic characteristics of MPXV isolated in Emilia-Romagna

Eleven MPXV strains, isolated, between May and September 2022 in Emilia Romagna, from 10 males and 1 female patients, were submitted for genomic characterization. All the MPXV samples were preliminarily tested through real-time PCR and confirmed as Clade IIb¹³. To explore the molecular evolution and genomic characteristics of MPXV associated with the Emilia-Romagna outbreak, isolates underwent whole-genome sequencing. An initial phylogenetic tree on Nextclade (https://clades.nextstrain.org) confirmed that all the analyzed isolates clustered with clade IIb. By integrating the 11 MPXV sequences obtained from Emilia-Romagna with a curated dataset of 267 high-quality Clade II MPXV genomes (Supplementary Table S1), phylogenetic analyses revealed that all strains isolated from male patients grouped within lineage B, being distributed across different sub-lineages: B.1 (n = 5), B.1.3 (n = 3), and B.1.12 (n = 2) (Fig. 1). Specifically, the sequences hMPXV/P1/2022, hMPXV/P4/2022, hMPXV/P9/2022, hMPXV/P10/2022, and hMPXV/P15/2022 clustered within lineage B.1 with strong bootstrap support (> 80). The sequences hMPXV/P3/2022, hMPXV/P5/2022, and hMPXV/P6/2022 were grouped within lineage B.1.3, supported by high bootstrap values (> 90). While hMPXV/P7/2022 and hMPXV/P13/2022 clustered within lineage B.1.12, although with lower bootstrap support (< 50). In contrast, the isolate hMPXV/P12/2022, collected from the only female patient, clustered within lineage A, sub-lineage A.2.3. with maximal bootstrap support (100). Epidemiological data indicated that this case was linked to travelling from Ghana to Italy (Table 1).

These results were further confirmed by the phylogenetic analysis performed using the Nextclade tool.

Phylogenetic Analysis of 11 MPXV Genomes from Emilia-Romagna. The Maximum Likelihood tree was built based on the whole-genome alignment of the 11 Italian MPXV sequences obtained in this study, alongside 267 high-quality clade IIb MPXV genomes retrieved from NCBI GenBank (as of October 4, 2024). The tree was rooted using three clade IIa sequences, with outgroup sequences highlighted in light brown. Tree labels are color-coded by lineage. Bold red labels in the tree indicate the 11 sequences from this study emphasized in bold red. Black dots indicate bootstrap support, with only values > 70 shown, while nodes with bootstrap values < 20 were collapsed. The scale bar represents the number of nucleotide substitutions per site. The tree was inferred using IQ-TREE v.2.1.4, with automated model selection (best fit model was HKY + F + R4 according to Bayesian Information Criterion) and 10,000 ultrafast bootstrap replicates. The middle annotation ring indicates the isolation location of each MPXV genome, while the outer annotation ring indicates the corresponding isolation date.

Sequences from lineage A and its related sub-lineages were retrieved and analysed in detail. As confirmed by sequences available in both NCBI GenBank and GISAID, most lineage A sequences have been isolated from the African continent, particularly Nigeria. In contrast, the few A.1.1 sequences, available in public databases, have been identified in a more diverse range of geographical regions, including the USA, Italy, and Thailand (Fig. 2).

Geographic distribution of MPXV A sub-lineages. The map illustrates the global distribution of MPXV A sub-lineages, highlighting their widespread presence, particularly in the southern provinces of Nigeria. The map was generated using QGIS v.3.34.14 (http://www.qgis.org). In cases where detailed location data was unavailable, points were placed at the centroid of the respective countries where the virus was detected. The sequences used for this analysis include all lineage A sequences downloaded from NCBI as of October 4th.

Temporal increase of SNPs and APOBEC3 associated mutations in MPXV lineage A

APOBEC analysis confirmed that early MPXV samples, isolated between 2017 and 2020, generally showed lower APOBEC mutation counts, while those isolated between 2021 and 2024 show an increased range of mutation counts, including higher values in the 15 to 35 range (Fig. 3A). As expected, A.2 and A.3 sub-lineages generally exhibit higher APOBEC-style mutations, mostly covering a range of 16 to 35 counts, against the range covered by lineage A and A.1 which is 5 to 15 (Fig. 3B).

APOBEC-associated mutations in MPXV sub-lineages over time. (A) Yearly distribution of APOBEC-associated mutations per sample across different sub-lineages. Circles indicate samples retrieved from NCBI (accessed on October 4, 2024), whereas the cross denotes the sample sequenced in this study. (B) Number of APOBEC-associated mutations in MPXV clade IIb A sub-lineages.

In the analyzed A.2 lineage sequences, we found on average 57.3 SNPs (σ = 24.1) per sample. The high standard deviation can be explained by the presence of two samples (GenBank Accession Number: PP853014 and PP853015) with high numbers of SNPs, 110 and 94 respectively. In the A.2.1 sub-lineage (seven samples: three from South Korea, two from the USA, one from Egypt and one from the UK) we found on average 44.4 SNPs (σ = 5.9) per sample. In the A.2.2 sub-lineage, we analyzed fourteen samples (one from the UK and thirteen from Nigeria), with an average mutation number of 72.7 (σ = 23.4). In the A.2.3 sub-lineage, 47 samples, one from Australia, one from the USA, one from Italy and 44 from Nigeria showed an average number of mutations of 56.3 (σ = 15.7).

Analysis of SNPs and APOBEC-associated mutations in sub-lineage A.2.3

To gain deeper insight into the evolution pathway of the hMPXV/P12/2022, we conducted a detailed examination of SNPs and APOBEC-associated mutations within the whole sub-lineage A.2.3. This analysis identified 36 unique polymorphisms across the genome, including 3 A/G, 15 G/A, and 18 C/T substitutions (Fig. 4A). Of the observed G-to-A transitions, 54.4% occurred in an APOBEC context, of these 90.6% were specifically 5’ GA-to-AA while the remaining 9.4% were GG-to-AG. On the other hand, of the observed C-to-T transitions, 56.9% occurred in an APOBEC context, meaning a TC-to-TT mutations. The prevalent 5’ APOBEC mutation is TC-to-TT in strains of lineage A, representing 51.5% of all the APOBEC-style mutations, followed by the GA-to-AA mutation (44.0%). Considering the subset of strains of lineage A.2.3, the TC-to-TT mutation is also in this case the most present, being 50.8% of the APOBEC mutations in this subset, followed by 44.2% of GA-to-AA mutations (Fig. 4B).

Genomic variability and APOBEC mutations in the A.2.3 sub-lineage. (A) Histogram representing the detected variants in the A.2.3 sub-lineage in our dataset. Variants are expressed with the genomic position, the expected nucleotide in the reference genome for lineage A and the alternate nucleotide detected in our samples. The x-axis represents the genomic position and the corresponding SNP type. The y-axis expresses the count both in terms of absolute occurrences (left) and in terms of percentage with respect to this experimental dataset (right). (B) Histogram representing the APOBEC Mutation counts per sample categorized by mutation type in the A.2.3 sub-lineage. The x-axis displays the sequences included in our dataset. The y-axis represents the absolute number of APOBEC-style mutations.

Signature mutations and gene disruption in A.2.3 sub-lineage

Four specific transitions were found in genomes belonging to lineage A.2.3 at positions 57,284 (C/T), 146,223 (C/T), 170,317 (G/A), and 175,759 (C/T). These are associated with alterations to the protein primary sequence. They are in the OPG075 (Glutaredoxin-1, R72K), OPG170 (Chemokine binding protein, D186N), OPG198 and OPG205 (Ser/thr kinase, D94N) and (Ankyrin repeat protein, A35V). Overall, hMPXV/12/2022 showed a total of 47 single SNPs relative to the reference MPXV-M5312_HM12_Rivers (GenBank Accession Number: NC_063383) causing alterations in the primary protein sequences ¹⁴. Of these, 3 were identified as strain specific: OPG023 leading to gene disruption at position 12,207 (G/A), OPG094 at position 74,420 (C/T) and OPG143 at position 122,712 (C/T) causing amino acid changes A214V and D202N respectively in the Entry/fusion complex component similar to VACV-Cop G9R and in a soluble myristilated protein similar to VACV-Cop A16L,

Another mutation was identified at OPG188 at position 161,953 (C/A) leading to an amino acid change (P314T). This nonsynonymous mutation appeared to be shared with 2 other strains of lineage A.2.3 isolated from Nigeria in 2022 and 2023 (GISAID: EPI ISL 19256279 and 19256278) and lead to a mutation in a schlafen-like protein similar to VACV-Cop B2R.

Concerning the MPXV replication complex, all the identified B sub-lineages of this study showed nonsynonymous mutations in OPG071 and OPG093 encoding respectively similar to the VACV-Cop F8L (RC catalytic subunit) and G9R (a processivity factor). In particular, a change at position 108 (F/L) of the encoded F8L and at position 30 (S/L) and 88 (D/N) of the G9R were observed. These mutations are ~ 100% prevalent in the strains that circulated during the 2022 pandemic. These mutations were not observed in the ancestor Nigeria 1971 (GenBank Accession Number: KJ642617.1) and lineage A (sub-lineages A.1, A.2, A.2.1, A.2.2, A.2.3) which conserves amino acid residues L at position 108 of the F8L and S and D at positions 30 and 88 of the G9R as identified in strains isolated by 2018. Sub-lineage A.1.1 seems to represent an exception, within the A lineage, as it resembles lineages B with F at position 108 of the F8L and L at position 30 of the G9R, while at position 88 of this protein A.1.1 shows the residue D as all the viruses of A lineage.

Two non-sense mutations were found in hMPXV/P12/2022 one of which was previously identified in other strains classified as A.2 sub-lineages. This nonsense mutation is in OPG176 encoding a IL-1/TLR signalling inhibitor similar to VACV-Cop A46R, that was reported to be compatible with APOBEC activity according to Ndodo et al.¹⁵. As previously described, this gene disruption is only found in lineage A.2. In fact, we were able to detect the mutation only in the strain A.2.3 hMPXV/P12/2022, but not in the strains belonging to B sub-lineages confirming that this genetic disruption is the consequence of hallmark adaptation event of the A.2 lineage¹⁶. This gene disruption leads to the translation of a polypeptide of 21 amino acids, instead of a whole-length protein of 240 residues. This variant was reported to be actively circulating among human hosts in West Africa and outside¹⁵.

A further gene disruption was found in the hMPXV/P12/2022 in the OPG023 encoding one of the MPXV Ankyrin repeat proteins similar to VACV-Cop D7L. The nucleotide G/A substitution at position 12,207 causes the appearance of a stop codon with subsequent protein truncation at amino acid position 232 (Q/*), which appeared to be a very rare variant, only found in another A.2.3 lineage strain detected in Nigeria in 2022 (GISAID: EPI ISL19256202). To characterize the previously unreported effects of the gene disruption in OPG023, we queried the reference database for protein search in UniProt¹⁷ (https://www.uniprot.org/). As expected, the 232 amino-acid residues truncated protein of hMPXV/P12/2022 (Fig. 5A) showed to be 100% identical to the Ankyrin repeat-containing protein PG023 of MPXV-M5312_HM12_Rivers (UniProt ID: A0A7H0DMZ9) which is overall 660 residues long. A 94% identity was identified with the D6L protein of the Variola virus isolate Human/India/Ind3/1967 (UniProt ID: Q07045). The protein superimposition with the VARV 452 residue long protein, modeled with the homology-based platform SWISS-MODEL (Fig. 5B) showed a high structural identity in the conserved ANK repeats, characterized by conserved domain of approximately 33 amino acids as originally identified in Ankyrin. These motifs have an L-shaped structure consisting of a beta-hairpin and two alpha-helices.

Structural modeling of hMPXV/P12/2022 and comparison with the D6L protein. (A) Structural AlfaFold3 model of hMPXV/P12/2022 The protein contains 12 Ankyrin repeats from the N to the C terminus. The computed average per-atom model confidence score is 84, indicating that residues are modeled on average with high quality (see “Methods”). (B) Structural model of D6L protein (Uniprot ID: Q07045). Model is computed by SWISS MODEL with an Average Model Confidence (QMEANDisCo) is equal to 0.58 ± 0.05 (see “Methods”). The protein length is 452 residues. The structural alignment of the models of hMPXV/P12/2022 with the D6L protein. The two models superimpose with a root mean square deviation (RMSD) of 1.65 Å over 216 aligned residues (sequence identity is 94%).