Genetic recombination in fast-spreading coxsackievirus A6 variants: a potential role in evolution and pathogenicity

Abstract

Hand, foot, and mouth disease (HFMD) is a common global epidemic. From 2008 onwards, many HFMD outbreaks caused by coxsackievirus A6 (CV-A6) have been reported worldwide. Since 2013, with a dramatically increasing number of CV-A6-related HFMD cases, CV-A6 has become the predominant HFMD pathogen in mainland China. Phylogenetic analysis based on the VP1 capsid gene revealed that subtype D3 dominated the CV-A6 outbreaks. Here, we performed a large-scale (near) full-length genetic analysis of global and Chinese CV-A6 variants, including 158 newly sequenced samples collected extensively in mainland China between 2010 and 2018. During the global transmission of subtype D3 of CV-A6, the noncapsid gene continued recombining, giving rise to a series of viable recombinant hybrids designated evolutionary lineages, and each lineage displayed internal consistency in both genetic and epidemiological features. The emergence of lineage -A since 2005 has triggered CV-A6 outbreaks worldwide, with a rate of evolution estimated at 4.17 × 10^-3 substitutions site^{- 1} year^-1 based on a large number of monophyletic open reading frame sequences, and created a series of lineages chronologically through varied noncapsid recombination events. In mainland China, lineage -A has generated another two novel widespread lineages (-J and -L) through recombination within the enterovirus A gene pool, with robust estimates of occurrence time. Lineage -A, -J, and -L infections presented dissimilar clinical manifestations, indicating that the conservation of the CV-A6 capsid gene resulted in high transmissibility, but the lineage-specific noncapsid gene might influence pathogenicity. Potentially important amino acid substitutions were further predicted among CV-A6 variants. The evolutionary phenomenon of noncapsid polymorphism within the same subtype observed in CV-A6 was uncommon in other leading HFMD pathogens; such frequent recombination happened in fast-spreading CV-A6, indicating that the recovery of deleterious genomes may still be ongoing within CV-A6 quasispecies. CV-A6-related HFMD outbreaks have caused a significant public health burden and pose a great threat to children's health; therefore, further surveillance is greatly needed to understand the full genetic diversity of CV-A6 in mainland China.

Keywords: Coxsackievirus A6; evolutionary dynamics; foot and mouth disease; genetic recombination; hand; pathogenicity; phylogeny; whole-genome analysis.

Figure 1.

Maximum-likelihood phylogenetic trees of 142 global CV-A6 variants constructed based on the (A) P1 capsid region and (B) noncapsid region; the lineages are differentiated by distinct colors, which are indicated in the bottom left, and each subtype is indicated on the right side of the P1 tree. The branches in each phylogeny are colored according to lineage. (C) Intralineage pairwise similarity comparison based on ORF sequences computed between each sequence and the group mean using sliding window nucleotide similarity analysis with a 200-nt window moving in 20-nt steps; similarity was not calculated for lineages containing only one sequence.

Figure 2.

The MCC phylogenetic tree generated using the MCMC method based on the entire P1 sequences of 142 global CV-A6 variants and colored according to different lineages. The scale bar represents time in years. The tree was node-labeled with inferred dates of lineage splits. Each subtype is shaded in light gray, except for the independent lineage –G. Recombination breakpoints based on the ORFs of different lineages were detected within (B) subtype D1, (C) subtype D2, and (D) subtype D3 using the original lineage of each subtype as the query group with a 200-nt window moving in 20-nt steps.

Figure 3.

(A) Time-scaled phylogenetic tree generated using the MCMC method for 336 complete CV-A6 P1 sequences from mainland China. The seven lineages from mainland China are indicated by different branch colors. The scale bar represents time in years. The tree was node-labeled with inferred dates of lineage splits. (B) The geographical distribution of seven lineages in mainland China. The proportion and number of cases caused by each lineage are indicated by the pie chart located in each geographic region on the map. The geographic map of China was taken from Highcharts (grant number: 0321912045738052), with regions denoted by different colors. Similarity plots were constructed and bootscan analyses were performed for parental recombination detection among screened EV-A strains and (C) lineage –K2, (D) lineage –J, and (E) lineage –L. Using lineage –A as the comparative group, prototype strains of CV-A6 and CV-A2 were used as outgroups. (F) A Gaussian Markov random field (GMRF) skyride plot of the P1 region of Chinese CV-A6, reflecting the relative genetic diversity from 2010 to 2018. The x-axis is the time scale (years), and the y-axis is the effective population size. The solid line indicates the median estimates, and blue shading indicates the 95 per cent HPD.

Figure 4.

Comparisons of pairwise distance in P1- and 3D^pol-region sequences within seven lineages; between lineage –A and lineages –J, L, and –K2 (from subtype D3) and between lineage –A and lineages –C,–D, and –K1 (from subtype D2). Axes depict nucleotide divergence between the two genomic regions. Each lineage is indicated by a differently colored dot in the scatterplot.

Figure 5.

MCC phylogenetic trees based on the ORF sequences of (A) lineage –A, (b) lineage –J, and (C) lineage –L from mainland China. The branches are colored according to the seven geographic regions of mainland China, which are indicated in the bottom left. The scale bar represents time in years. The tree was node-labeled with inferred dates of lineage splits.

Figure 6.

The aa polymorphisms in the P1 capsid region (A) between subtypes D2 and D3 and (B) among lineages –A, –J, and –L within subtype D3. The positions of aa polymorphisms detected when comparing subtype D3 with subtype D2 are marked in blue, those detected when comparing lineages –A, –J, and –L within subtype D3 are marked in crimson, and position 242 is marked in both colors. The positions of aa polymorphisms that are embedded in the assumed epitope loops are marked with a superscript asterisk. The web-based application WebLogo was used for generating aa sequences logos (Crooks et al. 2004) (http://weblogo.threeplusone.com/).