The association of recombination events in the founding and emergence of subgenogroup evolutionary lineages of human enterovirus 71

Abstract

Enterovirus 71 (EV71) is responsible for frequent large-scale outbreaks of hand, foot, and mouth disease worldwide and represent a major etiological agent of severe, sometimes fatal neurological disease. EV71 variants have been classified into three genogroups (GgA, GgB, and GgC), and the latter two are further subdivided into subgenogroups B1 to B5 and C1 to C5. To investigate the dual roles of recombination and evolution in the epidemiology and transmission of EV71 worldwide, we performed a large-scale genetic analysis of isolates (n = 308) collected from 19 countries worldwide over a 40-year period. A series of recombination events occurred over this period, which have been identified through incongruities in sequence grouping between the VP1 and 3Dpol regions. Eleven 3Dpol clades were identified, each specific to EV71 and associated with specific subgenogroups but interspersed phylogenetically with clades of coxsackievirus A16 and other EV species A serotypes. The likelihood of recombination increased with VP1 sequence divergence; mean half-lives for EV71 recombinant forms (RFs) of 6 and 9 years for GgB and GgC overlapped with those observed for the EV-B serotypes, echovirus 9 (E9), E30, and E11, respectively (1.3 to 9.8 years). Furthermore, within genogroups, sporadic recombination events occurred, such as the linkage of two B4 variants to RF-W instead of RF-A and of two C4 variants to RF-H. Intriguingly, recombination events occurred as a founding event of most subgenogroups immediately preceding their lineage expansion and global emergence. The possibility that recombination contributed to their subsequent spread through improved fitness requires further biological and immunological characterization.

Fig 1

Phylogeny of VP1 and 3Dpol regions of EV71 for study subjects and those previously determined (listed in Table 1; see also Table S1 in the supplemental material). Clades were identified by bootstrap analysis (values of ≥70%) by maximum likelihood analysis as implemented using RAxML. The sizes of the triangles are proportional to the number of sequences within each clade (bootstrap resampling values are shown on branches). By branch rotation to maximize visual congruence of the two trees and the use of bootstrap values of ≥70% to define phylogenetic groupings, the minimum numbers of incongruent phylogeny relationships (depicted by red dotted lines) were determined (clades showing congruent branching orders labeled with blue dotted lines). The 3Dpol region includes available sequences from CVA16; these are labeled red, as their interspersed positions in 3Dpol are invariably incongruent with their outgroup position in VP1 (not included in the tree on the left).

Fig 2

Relationships between sequence divergence (MCL pairwise distances) in the VP1 region (x axis) and in 3Dpol (y axis) among (A) genogroup B (A) and genogroup C (B) EV71 sequences. In the examples shown, sets of pairwise distances in both regions between GgB4 and C1 to other variants within GgB and GgC are depicted. Note that subgenogroups C2 and C4 both contained within them single sporadic recombinants; these account for the additional groupings of data points encircled by the dotted line.

Fig 3

Association between VP1 sequence divergence (the maximum value is shown for each bar on the x axis) and the proportion of recombinant comparisons (i.e., belonging to different 3Dpol clades) for GgB and GgC.

Fig 4

MCMC tree of the VP1 sequences of GgB (A) and GgC (B) from the Asian Pacific region visualized using FigTree and plotted on a temporal y-axis scale using their sampling dates. Branches are color coded (see the key in each panel) according to the recombination group of individual sequences and their reconstructed ancestors.