Genomic Epidemiology and Phylodynamic Analysis of Enterovirus A71 Reveal Its Transmission Dynamics in Asia

Abstract

Enterovirus A71 (EV-A71) is one of the main pathogens causing hand, foot, and mouth disease (HFMD) outbreaks in Asian children under 5 years of age. In severe cases, it can cause neurological complications and be life-threatening. In this study, 200 newly sequenced EV-A71 whole-genome sequences were combined with 772 EV-A71 sequences from GenBank for large-scale analysis to investigate global EV-A71 epidemiology, phylogeny, and Bayesian phylodynamic characteristics. Based on the phylogenetic analysis of the EV-A71 3D^pol region, six new evolutionary lineages (lineages B, J, K, O, P, and Q) were found in this study, and the number of evolutionary lineages was expanded from 11 to 17. Temporal dynamics and recombination breakpoint analyses based on genotype C revealed that recombination of nonstructural protein-coding regions, including 3D^pol, is an important reason for the emergence of new lineages. The EV-A71 epidemic in the Asia-Pacific region is complex, and phylogeographic analysis found that Vietnam played a key role in the spread of subgenotypes B5 and C4. The origin of EV-A71 subgenotype C4 in China is East China, which is closely related to the prevalence of subgenotype C4 in the south and throughout China. Selection pressure analysis revealed that, in addition to VP1 amino acid residues VP1-98 and VP1-145, which are associated with EV-A71 pathogenicity, amino acid residues VP1-184 and VP1-249 were also positively selected, and their functions still need to be determined by biology and immunology. This study aimed to provide a solid theoretical basis for EV-A71-related disease surveillance and prevention, antiviral research, and vaccine development through a comprehensive analysis. IMPORTANCE EV-A71 is one of the most important pathogens causing HFMD outbreaks; however, large-scale studies of EV-A71 genomic epidemiology are currently lacking. In this study, 200 new EV-A71 whole-genome sequences were determined. Combining these with 772 EV-A71 whole-genome sequences in the GenBank database, the evolutionary and transmission characteristics of global and Asian EV-A71 were analyzed. Six new evolutionary lineages were identified in this study. We also found that recombination in nonstructural protein-coding regions, including 3D^pol, is an important cause for the emergence of new lineages. The results provided a solid theoretical basis for EV-A71-related disease surveillance and prevention, antiviral research, and vaccine development.

Keywords: Bayesian phylodynamics; EV-A71; hand foot and mouth disease; phylogenetics; recombination.

FIG 1

Phylogenetic analysis of 259 EV-A71. (A) (Left) ML tree constructed from the VP1 region of 259 sequences with the corresponding genotype/subgenotype of each clade labeled; (right) ML tree constructed from the 3D^pol of 259 sequences. Different colors are used to distinguish the different lineages, and the two ML trees are colored according to the lineage to which the sequences belong. Dotted lines indicate subgenotypes and lineages with phylogenetic differences. (B) (Left) Heat map of the nucleotide similarity of the lineage E and H VP1 region sequences, where the two subgenotypes belonging to the same lineage have reduced similarity in the VP1 region, and different lineages have an even lower similarity; (right) heat map of nucleotide similarity of the C1, C2, and C4 3D^pol regions. Again, the similarity of different lineages belonging to the same subgenotype significantly decreases. (C) Pairwise distance comparison of VP1 and 3D^pol sequences. The area labeled “I” shows the comparison of intralineage sequences, and the rest compares different lineages with lineage I. The area labeled “II” compares lineages with lineage I of the same C4. The colors of different lineages are as in panel A. (D) Comparison of intralineage sequence similarity based on ORF region using sliding-window nucleotide similarity analysis, with 200 nucleotide windows moving in steps of 20 nucleotides. A lineage containing only 1 sequence was not analyzed.

FIG 2

Temporal dynamics analysis of global EV-A71. (A) A maximum clade credibility (MCC) tree was constructed based on 256 EV-A71 VP1 sequences, with the branches colored according to lineage and the tip colored according to region. (B to D) Temporal dynamics analysis of subgenotypes (left) and recombination breakpoint analysis (right) of subgenotypes C1, C2, and C4, respectively. The colors of the MCC tree are consistent with those in panel, and the query groups are lineage E, lineage E, and lineage I, respectively.

FIG 3

Phylogeographic analysis of B5 in Asia. (A) MCC tree constructed based on 85 B5 VP1 sequences. Branch colors indicate inferred location states, diamonds indicate the presence of B5 spread between regions, and the spread direction and time range are marked at the corresponding positions. (B) B5 spatial diffusion pathways in Taiwan (China), Vietnam, and Thailand. Only migration pathways with a mean indicator of >0.5 are shown. Solid black arrows indicate migration pathways with decisive support (BF > 1,000), dashed black arrows indicate migration pathways with very strong support (150 < BF < 1,000), solid gray arrows indicate migration pathways with strong support (20 < BF < 150), and dashed gray arrows indicate migration pathways with support (3 < BF < 20). All migration pathways in this paper follow this criterion. (C) Histogram of the total number of location state transitions inferred from the three regions.

FIG 4

Phylogeographic analysis of C4 in Asia. (A) MCC tree constructed based on 170 C4 VP1 sequences. Branch colors indicate inferred location states. Lineages N and O are highlighted by dots on the tip. Diamonds indicate the existence of C4 spread between regions, with the direction and time range of spread noted close to them. (B) Spatial diffusion pathways of C4 in China, Vietnam, and Cambodia, showing only migration pathways with BF values of >3 and indicators of >0.5. (C) Histogram of the total number of location state transitions inferred from the three regions.

FIG 5

Phylogeographic analysis of C4 in China. (A) Histogram of root state posterior probabilities in eight regions of China. Different colors represent different regions. The horizontal coordinate shows posterior probabilities, and the vertical coordinate shows regional distributions. (B) MCC tree constructed based on 234 VP1 sequences; colors of branches indicate inferred location states. Solid shapes represent the propagation of C4 between different regions. (C) Spatial diffusion pathways of C4 in eight regions in China. Only migration pathways with BF values of >3 and indicators of >0.5 are shown. (D) Histogram of the total number of inferred location state transitions for the eight regions.

FIG 6

Selective pressure analysis of EV-A71. The black line represents amino acid sequences translated by the EV-A71 VP1 region, and data sets from different regions and times are indicated by different colors. The presence of markers directly above and below the amino acid site indicated that the site had a positive selection. Squares indicate positive selection sites supported only by SLAC, circles indicate positive selection sites supported only by MEME, and stars indicate support by both methods. The markers directly above and below sites indicate positive selection sites found in data sets from different regions and years, respectively.