Transition of D3c branch and novel recombination events contribute to the diversity of Coxsackievirus A6 in Beijing, China, from 2019 to 2023

Abstract

Coxsackievirus A6 (CVA6) is a major pathogen responsible for numerous outbreaks of hand, foot, and mouth disease (HFMD) worldwide. This study investigates the molecular evolution and recombination of CVA6 in Beijing, China. Full-length sequences of 54 CVA6 from Beijing (2019-2023) were obtained through metagenomic next-generation sequencing and Sanger sequencing. These sequences were compared with representative sequences from GenBank to analyse their phylogenetic characteristics, recombination diversity, and evolutionary dynamics. The 54 CVA6 strains co-circulated with those from multiple provinces in China, as well as from South Korea and Japan. Phylogenetic analysis revealed a novel D3c branch, with the VP1 T283A amino acid mutation identified as a key change in its formation. One sequence belonged to the D3a branch, while 53 sequences belonged to the D3c branch. Recombination analysis identified RF-A (46, 85.1%) and three novel recombinant forms (RFs): RF-Z (1, 1.9%), RF-AA (1, 1.9%), and RF-AB (6, 11.1%). Bayesian phylogenetic analysis estimated that the most recent common ancestor of D3c emerged in August 2013 (95% highest probability density (HPD): May 2012 to September 2014), with recombination events occurring in RF-Z (2017-2019), RF-AA (2019-2023), and RF-AB (2021-2023). In conclusion, we revealed a globally circulating CVA6 D3c branch and identified three novel RFs, providing valuable insights for the intervention and control of HFMD.

Keywords: Coxsackievirus A6 (CVA6); hand, foot, and mouth disease (HFMD); metagenomic next-generation sequencing; novel D3c branch; recombinant forms (RFs).

Figure 1

The D3c branch is the predominant evolutionary branch. (A) Phylogenetic analysis of VP1 sequences (nucleotide positions 2441–3355, 915 nt, relative to AY421764) using 172 CVA6 sequences from Sequence Library 1. ML tree was constructed with 1000 bootstrap replicates to assess the robustness of the groupings, with bootstrap values ≥80% indicated. (B) Entropy values at each position of the 4736 VP1 amino acid sequences from Sequence Library 2. (C) Composition of VP1 residues at positions 30, 137, and 283 in the D3a, D3b, and D3c branches. (D) Atomic model of the CVA6 particle capsid, obtained from the Protein Data Bank (PDB ID: 5YHQ), visualized using PyMOL (version 2.6). VP0 is shown in raspberry, VP1 in sky blue, VP3 in tangerine orange, and VP1-283A in green. The C-terminus of VP1 is labelled in blue. VP1-283A forms hydrogen bonds with VP3-95Q and VP3-66R. (E) Yearly distribution of the CVA6 D3 sub-genotype (D3a, D3b, and D3c from Sequence Library 2) worldwide from 2008 to 2024 (top). Yearly distribution of the CVA6 D3 sub-genotype (D3a, D3b, and D3c from Sequence Library 3) in Beijing from 2013 to 2023 (bottom).

Figure 2

Three novel recombination forms in eight CVA6 sequences. (A) Phylogenetic comparison of partial 3Dpol sequences (nucleotide positions 6073–6884, 812 nt, relative to AY421764) using 54 sequences and representative CVA6 sequences obtained from GenBank. ML tree was constructed with 1000 bootstrap replicates to assess the robustness of the groupings, with bootstrap values ≥80% shown. Recombination events from (B) CVA6/C-Tan/BJ/202301, (C) CVA6/C-Tan/BJ/201907, and (D) CVA6/C-Tan/BJ/202308 sequences were identified through similarity plots and bootscanning analyses, using a sliding window of 200 nt with 20 nt steps.

Figure 3

Genomic map of recombination events in eight CVA6 sequences predicted by RDP4. The black band represents the full-length genome of the CVA6, with numbers above indicating the start and end positions of the breakpoints (numbers in brackets represent the breakpoint positions aligned with AY421764). The grey bands represent the recombined genomic regions, with the text below indicating the major and minor parents of the predicted recombination event.

Figure 4

Parental strains identification for RF-AA in the 6156–7321 nt region. (A) Analysis of the nucleotide and amino acid sequence similarity of CVA6/C-Tan/BJ/202308 compared to the complete genomic sequences of CVA2, CVA4, CVA6, and CVA8 available in GenBank. (B) Phylogenetic trees of RF-A and the 10 most similar sequences, estimated from the 6155–7297 nt region (top) and synonymous sites only (bottom). Clade support values obtained from 1000 bootstrap replicates are shown.

Figure 5

Analysis of mutation profiles in the RF-Z, RF-AA, and RF-AB sequences. The numbers above each site correspond to the positions of the amino acid with KM114057 (the earliest reference strain of sub-genotype D3). Recombinant regions associated with the RF-Z and RF-AB sequences are highlighted by black boxes. The 6155–7297 nt region in the RF-AA sequence is highlighted by an orange box. The VP1 surface loops are shaded in grey.

Figure 6

Global evolutionary dynamics of CVA6. (A) Bayesian skyline plot of the CVA6 VP1 region sequence, reflecting the relative genetic diversity of CVA6 from 2010 to 2024. The X-axis represents the time scale (in year), and the Y-axis shows the effective population size. The solid line is the estimated mean, and the blue shadow is the 95% highest posterior density interval. (B) The MCC phylogenetic tree generated using the MCMC method based on complete VP1 sequences of 172 CVA6. The colour of the branches represents the geographical location of sequences.